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Chapter  1 


Introduction  and  Background 

Non-smooth  systems  can  be  generally  defined  as  systems  where  the  vector  field 
representation  of  the  system  is  non-smooth.  A  typical  state-space  representation  of 
a  system’s  dynamics  is  shown  in  Eq.  (1.1),  where  a;  is  a  vector  of  the  system  states, 
t  is  time,  and  u  represents  a  vector  of  the  system  inputs.  By  using  the  state-space 
representation,  here,  a  non-smooth  system  is  defined  as  a  system  where  the  function 
/  is  not  Cl  continuous  with  respect  to  time  and  the  states. 

x  =  /(*,x,u(t))  (1.1) 

Non-smoothness  can  occur  due  to  non-smooth  forces,  system  geometry,  or 
non-smooth  inputs.  When  the  discontinuities  that  cause  the  non-smoothness  are 
finite  and  well  defined,  it  is  often  helpful  to  treat  the  non-smooth  system  as  a 
switched  combination  of  smooth  systems.  In  most  systems  of  interest  (in  terms 
of  the  importance  of  including  the  non-smoothness),  the  dynamics  of  the  switched 
systems  are  dramatically  different.  Typical  examples  of  non-smoothness  are  systems 
with  friction  or  impacts. 

In  general,  vehicle  systems  are  systems  wherein  a  vehicle  moves  in  space.  Often 
times,  the  forces  acting  upon  the  vehicle  are  best  defined  in  a  moving  reference  frame 
that  follows  the  vehicle.  This  is  because  in  these  systems,  many  of  the  forces  relate 
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Figure  1.1:  Supercavitating  vehicle. 

to  relative  positions  and  angles  with  respect  to  the  vehicle.  An  example  would  be 
tire  forces  for  a  wheeled  terrestrial  vehicle  that  depend  heavily  on  slip  angles,  and 
these  are  easiest  defined  with  respect  to  the  vehicle  orientation.  Another  example 
is  from  flight  dynamics,  wherein  the  lift  and  drag  forces  generated  from  the  control 
surfaces  depend  on  the  angle  of  attack. 

1.1  Supercavitating  Vehicles 

The  physical  problem  of  interest  for  this  work  is  supercavitating  underwa¬ 
ter  vehicles.  Cavitation  is  generally  described  as  vapor  bubble  formation  caused  by 
pressure  drops  associated  with  a  body  moving  in  a  fluid.  Cavitation  is  generally  con¬ 
sidered  as  being  a  detrimental  and  corrosive  process,  since  the  effects  of  collapsing 
vapor  bubbles  causes  damage  to  the  body  in  the  fluid.  Collapsing  bubble  cavitation 


2 


is  known  as  initial  cavitation,  and  this  type  of  cavitation  is  widely  associated  with 
pumps  and  propellers  [5].  However,  partial  cavitation  and  supercavitation,  relate 
to  sustained  cavities.  Partial  cavitation  relates  to  a  sustained  cavity  that  encom¬ 
passes  some  part  of  the  body  and  supercavitation  relates  to  a  sustained  cavity  that 
encompasses  the  entire  body  [40].  This  leads  to  supercavitating  vehicles,  which  are 
high-speed  underwater  vehicles  wherein  a  gas  cavity  surrounds  the  entire  body  of 
the  vehicle.  For  these  vehicles,  the  cavity  is  formed  by  vapor  due  to  cavitation  at  the 
nose  (created  by  a  blunt  trailing  edged  cavitator,  Figure  1.1)  and  this  cavity  may 
be  aided  through  forced  ventilation.  Supercavitating  vehicles  have  the  advantage  of 
reduced  drag  (due  to  the  reduction  in  wetted  surface  area),  which  allows  for  high 
speeds.  The  high-speed  capabilities  are  particularly  advantageous  for  weapon  or 
counter-weapon  systems,  although  other  applications  such  as  high-speed  underwa¬ 
ter  transport  have  been  proposed.  Supercavitating  vehicles  face  stability  challenges 
due  to  the  nonlinear,  non-smooth  planing  forces  that  arise  from  interaction  between 
the  body  and  the  cavity  wall. 

If  the  planing  force  function  with  respect  to  immersion  is  smooth  (aside  from 
at  zero  immersion),  then  the  supercavitating  vehicle  system  can  be  considered  as 
a  piecewise  smooth  system.  When  considering  motions  in  the  vertical  plane,  in 
a  general  sense,  the  system  dynamics  can  be  divided  into  three  regions  with  the 
dynamics  in  each  region  individually  characterized  by  smooth  dynamics.  One  while 
planing  on  the  top  surface,  one  while  planing  on  the  bottom  surface,  and  one  where 
the  vehicle  is  completely  inside  the  cavity  and  is  touching  neither  surface.  These 
regions  are  illustrated  in  Figure  1.2.  For  the  most  part,  the  cavity  will  form  in 
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a)  Planing  on  top  surface 


c)  No  planing 


Flow 

Direction 


Figure  1.2:  Regions  of  smooth  dynamics. 

the  direction  of  the  flow  with  respect  to  the  cavitator  (the  cavity  shift  effects  due 
to  cavitator  angle  of  attack  can  also  play  a  role  and  this  is  discussed  later  in  this 
work).  The  switching  points  between  regions  are  determined  by  the  point  where  the 
body  begins  to  contact  the  cavity.  For  a  constant  cavity  shape,  the  switching  point 
is  defined  by  the  angle  made  by  the  body  with  respect  to  the  flow.  Similarly,  when 
assuming  a  constant  cavity  shape,  the  planing  force  can  also  be  written  as  a  function 
of  the  angle  of  the  body  with  respect  to  the  flow,  since  this  angle  also  determines  the 
immersion  area.  Constant  cavity  approximations  allow  for  simpler  characterization 
of  the  system  dynamics  and  they  are  often  utilized.  However,  within  this  work,  non¬ 
constant  cavity  shapes  are  considered.  This  inclusion  creates  not  only  non-constant 
switching  boundaries,  but  also  non-constant  switched  dynamics,  since  the  immersed 
area  (and  hence  planing  force)  depends  on  the  cavity  shape. 

Supercavitating  vehicle  systems  differ  from  the  traditional  impacting  non- 


4 


smooth  systems.  The  use  of  non-constant  cavity  shapes  creates  moving  switching 
points  and  non-constant  switched  dynamics.  In  the  supercavitating  vehicle  system, 
the  non-smoothness  is  not  directly  position  dependent  as  with  impacting  systems 
(instead  related  to  the  attitude  of  the  vehicle  with  respect  to  the  cavity).  Addition¬ 
ally  the  planing  is  “soft”  in  that  movements  well  into  the  cavity  surface  are  allowable. 
The  planing  force  formulations  can  also  be  highly  complex  and  nonlinear,  and  these 
formulations  do  not  lend  themselves  to  a  direct  analytical  treatment. 

1.2  Literature  Review 

1.2.1  Dynamics  and  Control  of  Non-smooth  Systems 

Much  work  has  been  conducted  in  the  area  of  non-smooth  system  dynam¬ 
ics,  particularly,  related  to  impact  and  friction  problems.  The  following  is  a  brief 
summary  of  this  work. 

Tao  and  Lewis  [41]  applied  adaptive  control  techniques  for  systems  with  non¬ 
smooth  nonlinearities  such  as  backlash,  dead- zone,  component  failure,  friction,  hys¬ 
teresis,  and  saturation  and  time-delay  effects.  Adaptive  control  is  particularly  at¬ 
tractive  for  systems  where  non-smoothness  is  unknown.  Neural  networks  and  other 
adaptive  techniques  are  used  for  fault  detection,  feedback  control,  friction  compen¬ 
sation,  and  control  of  linear  time-delay  systems. 

Queiroz,  Malisoff,  and  Wolensk  [38]  presented  a  collection  of  research  results 
dealing  with  non-smoothness  in  optimal  control  problems.  Non-smoothness  can 
arise  from  the  system  itself,  or  the  min  operations  carried  out  for  the  optimization. 
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This  research  survey  dealt  with  developments  in  mathematical  control  theory,  opti¬ 
mization,  and  control  engineering  brought  on  by  advances  in  non-smooth  analysis. 
New  optimality  conditions  for  delay  differential  inclusions  and  effects  of  non-smooth 
analysis  on  Hamilton- Jacobi  and  Lyapunov  function  theory  have  been  presented. 

Li,  Soli,  and  Wen  [24]  look  at  switched  systems  that  are  defined  as  sets  of  con¬ 
tinuous  variable  systems  with  a  discrete  event  system  that  controls  the  switching. 
These  systems  can  be  characterized  as  having  non-smooth  behavior.  They  present 
a  cycle  analysis  method  and  some  results  about  conditions  for  Lyapunov  functions 
within  each  type  of  cycle  and  continuous  system.  A  method  to  provide  chaotic 
synchronization  for  chaos  based  encryption  systems  is  also  presented.  Fillipov  and 
Aizerman  theories  for  piecewise  continuous  systems  are  presented  in  reference  [3]. 
Numerical  schemes  for  convergence,  and  bifurcations  and  chaos  of  a  van  der  Pol- 
Duffing  oscillator  with  columb  friction  are  presented  along  with  several  other  system 
examples.  A  collection  of  lecture  notes  for  non-smooth  dynamical  systems  is  pre¬ 
sented  in  reference  [21].  Presented  are  Lyapunov  exponents  for  non-smooth  systems, 
Conley  index  theory,  KAM  theory,  and  Melnikov’s  method.  An  extension  of  Mel¬ 
nikov’s  method  is  presented  in  reference  [2],  Here,  the  method  is  extended  to  study 
higher  dimensional  systems  (greater  than  three).  They  are  able  to  predict  chaotic 
orbits  for  simple  frictional  models  (stick-slip  models),  and  attempt  to  apply  the  Mel¬ 
nikov  method  to  other  physical  systems  such  as  those  with  multi-body  dynamics. 
Lcine  and  Nijmeijer  use  convex  analysis  to  look  at  the  dynamics  and  bifurcations  of 
non-smooth  systems  [23].  This  combines  the  non-smooth  mechanics  approach  with 
non-smooth  dynamics  analysis. 
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Experimental  work  in  the  area  of  non-smooth  dynamics  has  also  been  con¬ 
ducted  in  studies  [33]  and  [36].  These  works  relate  to  bifurcations  in  systems  with 
impacting  dynamics.  Another  study  pertaining  to  cantilever  tip  impacts  is  presented 
in  reference  [6].  Drill-string  dynamics,  which  include  stick-slip  friction  interaction 
and  impacting,  is  discussed  in  reference  [25]. 

Much  of  the  work  in  the  area  of  mechanical  non-smooth  systems  relates  to 
friction  and  impact  type  dynamics.  Again  these  studies  differ  from  supercavitat- 
ing  systems  due  to  the  character  of  the  switching  boundaries,  the  relatively  “soft” 
switched  dynamics,  and  the  complexity  of  the  planing  forces. 

1.2.2  Trajectory  Planning  for  Vehicle  Systems 

Maneuvering  or  trajectory  planning  is  a  question  that  is  inherent  to  vehicle 
systems.  Since  motion  or  specific  maneuvers  for  these  systems  can  generally  be 
achieved  by  multiple  (or  potentially  infinite  number  of)  control  input  functions,  op¬ 
timal  control  methods  become  important  for  determining  “best”  inputs.  An  optimal 
control  approach  for  a  flexible  hull  swimming  vehicle  is  considered  in  reference  [4], 
Due  to  the  complexity  of  the  hydrodynamics  for  flexible  hull  vehicles,  an  analytic 
treatment  is  intractable.  Instead  a  numeric  genetic  evolutionary  algorithm  is  applied 
to  effectively  increase  swimming  performance.  An  optimal  approach  for  yaw  control 
on  a  four-wheel  vehicle  is  presented  in  study  [10].  Here,  parameters  of  a  particular 
control  form  are  optimized  to  enhance  yaw  performance.  A  hybrid  control  strategy 
for  an  autonomous  vehicle  is  considered  in  reference  [12],  Here,  feasible  trajecto- 
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ries  are  carried  out  by  lower  level  continuous  control  systems  where  a  higher  level 
control  layer  is  used  to  choose  between  discrete  feasible  trajectories  in  order  to  opti¬ 
mize  overall  performance.  This  type  of  optimization  can  be  applied  to  systems  with 
complex  dynamics  since  only  families  of  feasible  trajectories  need  to  be  considered. 
Motion  planning  is  also  discussed  with  the  option  of  including  obstacles  or  “no- 
fly”  areas.  Another  similar  hybrid  approach  is  used  in  reference  [11].  Again,  since 
only  a  discrete  set  of  feasible  trajectories  are  being  considered  by  the  upper  level 
control  scheme,  the  optimal  control  problem  can  be  solved  in  real  time.  Another 
computationally  efficient  optimal  control  problem  is  posed  in  reference  [17],  wherein 
trajectory  planning  is  considered  for  an  omni-directional  vehicle.  Here,  a  restricted 
set  of  allowable  controller  inputs  is  considered  to  reduce  optimization  complexity. 
Path  following  techniques  for  underwater  vehicles  are  presented  in  references  [22,  9]. 
In  these  studies,  the  vehicle  attempts  to  follow  a  prescribed  path. 

1.2.3  Supercavitating  Vehicle  Studies 

Dzielski  and  Kurdila  present  a  four-state  dynamic  model  of  a  supercavitating 
vehicle  [8].  A  simplified  cavity  model  and  planing  force  model  are  utilized.  With 
this  model,  unstable  and  limit-cycle  behavior  can  be  observed,  demonstrating  the 
complex  nature  of  the  vehicle  dynamics.  A  specific  linearizing  controller  was  also 
considered  in  which  the  nonlinearities  were  canceled  leaving  a  controllable  linear 
system.  This  approach  is  however  specific  to  the  model  parameters  and  assumes 
exact  knowledge  of  the  system  behavior.  Lin,  Balachandran,  and  Abed  [30]  looked 


at  the  use  of  a  switching  controller  for  the  same  system.  They  also  investigated 
the  presence  of  bifurcation  behavior  in  the  dynamic  system  and  presented  a  control 
method  to  delay  an  observed  Hopf  bifurcation  which  represents  the  onset  of  limit- 
cycle  (tail  slap)  behavior  [29].  In  reference  [26],  Lin  et  al.  present  results  on  absolute 
stability  for  sector  bounded  nonlinearities,  and  these  general  results  are  applied  to 
the  supercavitating  vehicle  system.  Kirschner,  Rosenthal  and  Uhlrnan  [20]  also 
looked  at  supercavitating  vehicle  dynamics.  In  their  study,  the  planing  force  is 
modeled  as  the  force  derived  from  a  nonlinear,  non-smooth  spring  and  damper 
combination.  In  this  work,  the  authors  also  accounted  for  the  time  delay  generated 
by  the  fact  that  the  cavity  radius  in  the  area  of  planing  is  actually  a  function  of 
previous  cavitator  position.  The  effects  of  the  time  delay  are  closely  examined  in 
references  [15]  and  [16].  In  these  studies,  control  efforts  that  work  to  stabilize  the 
delay-free  system  are  found  to  be  ineffective  in  stabilizing  the  delayed  system.  And 
for  certain  operating  conditions,  the  delay  has  a  destabilizing  effect.  Guidance  using 
a  tracking  method  is  considered  in  reference  [43].  Simple  tracking  maneuvers  are 
considered  by  using  two  outer-loop  control  strategies  consisting  of  a  pole  placement 
method  and  receding  horizon  control  scheme,  both  of  which  are  applied  to  a  feedback 
linearized  system.  Trajectory  optimization  for  supercavitating  vehicles  is  presented 
in  references  [39]  and  [1].  In  these  studies,  optimal  control  methods  are  used  to 
determine  the  maneuvering  characteristics  of  supercavitating  vehicles.  The  time- 
delay  effect  is  accounted  for  in  this  work,  along  with  fin  force  effects  generated 
by  the  distinct  individual  fin  immersion  depths  (clue  to  the  position  of  the  vehicle 
inside  of  the  cavity).  The  resulting  maneuvers  considered  in  this  work  do  not  involve 
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any  planing  of  the  vehicle.  In  the  previous  studies,  a  multitude  of  interesting  and 
complicated  dynamic  behavior  have  been  presented.  However,  the  cavity  and  planing 
force  models  are  often  simplified  and  this  is  one  area  of  research  intended  to  be 
pursued  in  this  work. 

There  has  also  been  work  in  the  area  of  cavity  modeling.  Initially  cavity  shapes 
were  used  to  design  body  shapes  along  constant  pressure  surfaces  [40].  There  have 
been  several  analytical  models  developed  for  predicting  cavity  shapes.  Logvinovich 
has  presented  a  detailed  analytical  discussion  of  cavity  shape  and  cavity  body  in¬ 
teractions  in  reference  [31].  An  iterative  potential  flow  model  for  supercavitation 
is  presented  in  studies  [42]  and  [19].  These  studies  help  generate  a  numeric  model 
for  cavity  shape  given  cavitator  and  flow  parameters.  Partial  cavity  dynamics  are 
significantly  more  difficult  to  model  since  the  body  protrusion  greatly  affects  the 
cavity  shape.  Several  schemes  are  outlined  in  reference  [40].  However,  all  of  these 
methods  for  approximating  partial  cavity  shapes  are  based  on  supercavitating  cavity 
models.  A  numeric  potential  flow  model  for  partial  cavity  shapes  is  presented  in 
reference  [44], 

Another  area  of  development  for  supercavitating  vehicles  is  with  the  planing 
forces  generated  by  the  body-cavity  interactions.  Analytical  planing  force  models 
are  presented  in  references  [32],  [37],  and  [14].  Since  planing  for  supercavitating 
vehicles  deals  with  non-flat  surfaces,  experimental  testing  is  complicated  and  data 
are  limited.  Some  model  validations  have  been  conducted  [7]. 


10 


1.3  Contributions 


The  first  aim  of  this  work  to  better  characterize  the  non-smoothness  in  the 
supercavitating  vehicle  system  and  explore  the  resulting  dynamics.  In  all  of  the 
previous  research,  cylindrical  cavities  without  cavitator  angle  of  attack  shift  effects 
are  assumed  in  order  to  simplify  planing  force  calculations.  In  this  work,  a  method 
to  allow  for  non-cylindrical  shifted  cavities  is  presented.  This  introduces  more  phys¬ 
ically  realistic  force  modeling,  which  however  dramatically  changes  the  nature  of 
the  non-smooth  forces.  This  creates  a  “soft”  non-smoothness  with  non-constant 
switching  boundaries  and  switched  dynamics.  Dynamics  and  bifurcations  of  system 
responses  are  explored  for  these  systems  as  well  as  stabilization  techniques. 

A  partially  cavitating  vehicle  dynamics  model  is  presented  in  an  effort  to  move 
towards  full  mission  modeling,  from  launch  (no  cavitation),  to  partial  cavitation 
(partially  wetted  vehicle),  to  full  supercavitation.  This  is  also  an  area  of  research 
not  accounted  for  in  any  of  previous  research  studies.  Without  the  planing  forces  at 
the  rear  of  the  vehicle,  it  is  found  that  a  more  active  fin  input  is  required  to  stabilize 
the  vehicle. 

Another  item  overlooked  in  much  of  the  previous  research  is  non-constant 
planing.  In  an  effort  to  reduce  the  complexity  of  the  planing  formulation,  steady 
planing  force  models  (not  accounting  for  vehicle  motions  into  or  out  of  the  cavity) 
are  often  utilized.  In  this  work,  a  model  that  accounts  for  these  “impacting”  forces 
is  presented.  Within  this  model  development,  a  full  description  of  the  origins  of  the 
immersion  depth  and  immersion  rate  terms  has  been  provided  for  use  in  the  planing 
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force  formulation.  This  is  an  aspect  that  has  been  overlooked  in  much  of  the  previous 
literature,  and  these  terms  have  been  previously  presented  in  a  somewhat  arbitrary 
manner.  Through  careful  representation  of  the  immersion  terms,  time-delay  effects 
can  also  be  properly  included  in  this  impacting  model.  The  time  delay  is  introduced 
by  the  fact  that  the  cavity  at  the  rear  of  the  vehicle  (where  planing  is  present), 
is  generated  by  a  previous  cavitator  position  and  orientation.  In  this  delayed  and 
impacting  model,  the  time  delay  is  found  to  be  stabilizing  for  certain  cavity  sizes. 
This  is  a  observation  that  is  contradictory  to  the  findings  of  previous  time-delay 
studies. 

All  of  these  modeling  efforts  are  meant  to  more  properly  characterize  the 
physical  non-smoothness  present  in  these  systems.  In  addition  to  the  dynamics  and 
stability  observations  made  for  these  systems,  maneuvering  is  also  considered  in  this 
dissertation  effort.  Maneuvers  are  presented  by  defining  specific  ending  conditions 
given  a  particular  initial  condition.  A  numeric  optimal  control  approach  is  utilized 
to  find  control  inputs  that  accomplish  the  maneuvers  in  the  quickest  possible  time. 
Unlike  other  optimal  control  studies  for  supercavitating  vehicles  [1,  39],  maneuvers 
in  this  work  are  considered  for  operation  regions  where  non-smooth  interactions  are 
prevalent  (such  as  speeds  where  the  cavity-body  clearances  are  tight).  The  resulting 
maneuvers  inherently  include  planing.  This  introduces  a  great  deal  of  complexity, 
especially,  since  the  instabilities  that  these  vehicle  systems  demonstrate  stem  from 
the  vehicle-cavity  interactions.  An  inner-loop  control  scheme  and  an  outer-loop  con¬ 
trol  scheme  are  proposed,  wherein  a  inner  feedback  loop  helps  stabilize  the  planing 
motions  of  the  vehicle,  while  the  outer-loop  guides  the  vehicle  towards  the  desired 
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end  condition.  This  control  configuration  works  well  with  the  inner-loop  rejecting 
the  fast  timescale  instabilities  not  addressable  by  the  coarse  outer-loop  control.  The 
non-smoothness  in  this  system  is  fairly  complex,  with  dramatically  different  dynam¬ 
ics  when  planing  (versus  not  planing),  and  small  windows  of  operation  (due  to  tight 
body-cavity  clearances)  which  leads  to  rapid  switching  between  different  dynam¬ 
ics.  With  such  complicated  dynamics,  simple  search  algorithms  that  are  coupled 
with  optimization  constraints  (determined  by  the  direct  integration  of  the  dynam¬ 
ics)  handled  as  penalties,  work  the  best  for  solving  the  optimal  control  problem. 
This  approach  could  be  extended  for  use  for  other  complicated  non-smooth  vehicle 
systems. 

1.4  Dissertation  Objectives  and  Organization 

The  overall  goal  of  the  work  is  to  achieve  a  better  understanding  of  the  dynam¬ 
ics  of  supercavitating  vehicles,  and  by  extension,  illustrate  a  method  to  analyze  the 
dynamics  of  other  non-smooth  systems.  Specific  objectives  include  the  following: 

1.  Examine  the  dynamics  of  supercavitating  vehicles  by  considering  cavity  models 
that  incorporate  realistic  physical  effects 

2.  Examine  the  dynamics  of  partially  cavitating  vehicles 

3.  Examine  the  dynamics  of  a  time  delayed,  non-constant  planing,  supercavitat¬ 
ing  vehicle  system 
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4.  Examine  the  maneuvering  capabilities  for  supercavitating  vehicles  and  pro¬ 
vide  a  framework  to  evaluate  maneuvering  for  other  complicated,  non-smooth 
vehicle  systems. 

The  rest  of  this  dissertation  is  organized  as  follows.  In  Chapter  two,  dynamics 
of  systems  using  shifted  cylindrical  cavities  are  discussed.  This  includes  a  description 
of  the  the  basic  dynamics  model  as  well  as  a  planing  force  formulation  that  is  used 
throughout  the  rest  of  the  work.  In  Chapter  three,  a  method  for  accounting  for 
non-cylindrical  cavities  is  introduced,  along  with  dynamics  results  for  these  systems 
that  can  be  compared  with  previous  research  hirelings.  Included  are  non-smooth 
bifurcations  and  methods  of  stabilizing  the  vehicle  motions  for  the  non-cylindrical 
cavity  system.  Additional  modeling  aspects  are  addressed  in  the  following  chapter. 
Here,  vehicle  dynamics  for  partial  cavities,  and  an  impacting  time-delayed  model  are 
presented.  Maneuvering  on  the  basis  of  an  optimal  control  approach  is  discussed  in 
Chapter  five.  Here,  maneuvers  are  considered  for  the  previously  discussed  systems, 
as  well  as  for  a  six  degree-of-freedom  system,  by  using  a  numeric  optimal  control 
approach  that  has  been  tailored  to  deal  with  the  non-smooth  vehicle  dynamics.  This 
is  followed  by  a  summary  of  the  contributions  and  a  discussion  of  potential  future 
directions.  Appendices  related  to  this  dissertation  are  included  at  the  end  to  provide 
additional  details  and  representative  algorithms  used  in  the  present  work 
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Chapter  2 


Vehicle  Dynamics  with  Cylindrical  Shifted  Cavities 

In  much  of  the  previous  vehicle  system  analysis,  the  cavity  is  assumed  to  be 
of  cylindrical  form  [1,  8,  15,  16,  29,  30,  39].  In  these  studies,  a  closed-form  solution 
for  the  cavity  radius  at  the  rear  of  the  vehicle  is  utilized  to  generate  an  approximate 
cylindrical  cavity  for  planing  force  calculations.  In  these  formulations,  the  influence 
of  the  cavitator  angle  of  attack  on  the  cavity  radius  are  ignored.  In  this  chapter,  a 
different  cavity  model  is  presented  that  allows  for  the  inclusion  of  cavitator  angle-of- 
attack  effects,  which  are  incorporated  as  a  shift  to  the  nominal  cavity  radius.  The 
cylindrical  assumption  is  maintained  in  order  to  highlight  the  cavity  shift  effects. 
A  more  accurate  shifted  non-cylindrical  cavity  planing  (utilizing  the  entire  cavity 
profile)  is  then  considered  in  the  following  chapter. 

To  help  illustrate  the  different  approaches  for  modeling  vehicle-cavity  inter¬ 
actions,  three  different  scenarios  are  illustrated  in  Figure  2.1.  In  part  a),  a  vehicle 
orientation  is  shown  with  a  typically  assumed  cylindrical  cavity.  The  cylinder  ra¬ 
dius  is  determined  by  an  estimate  of  the  nominal  cavity  radius  at  the  rear  of  the 
vehicle.  The  vehicle  body  is  depicted  to  be  entirely  within  the  cavity.  The  same 
body  orientation  with  a  shifted  cylindrical  cavity,  is  shown  in  part  b).  As  illustrated, 
the  cavitator  anglc-of-attack  can  shift  the  cavity  location  at  the  rear  of  the  vehi¬ 
cle  which  affects  when  planning  occurs  and  how  much  immersion  is  present  when 
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c)  Non-cylindrical  w/ Shift 


Figure  2.1:  Three  different  scenarios  in  modeling  body-cavity  interactions:  a)  vehicle 
within  a  cylindrical  cavity,  b)  vehicle  within  a  cylindrical  cavity  with  shifted  axis, 
and  c)  vehicle  within  a  non-cylindrical  cavity. 

planing.  A  shifted  full  non-cylindrical  cavity  is  shown  in  part  c).  The  immersion 
at  the  rear  of  the  vehicle  is  the  same  as  with  the  shifted  cylindrical  cavity  (since 
the  shifted  cylindrical  cavity  radius  and  position  is  determined  by  the  size  and  po¬ 
sition  of  the  non-cylindrical  cavity  at  the  rear).  However,  the  immersed  volume  is 
different,  which  changes  the  planing  forces  due  to  immersion. 

2.1  Basic  Dive-Plane  System  Model 

The  dive-plane  system  model  described  in  this  section  introduces  a  basic  vehi¬ 
cle  dynamics  representation  that  is  used  as  a  basis  for  the  subsequent  system  models 
presented  in  this  work.  The  basic  dive-plane  system  model  comes  from  references 
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X 


Figure  2.2:  Coordinate  system  definition  for  system  model. 

[29]  and  [30],  and  this  model  is  based  on  a  model  originally  presented  in  reference 
[8].  The  model  is  comprised  of  four  states  defined  in  a  non-inertial  reference  frame 
which  accounts  for  two-dimensional  ridged  body  movements  in  the  dive-plane.  The 
dive-plane  is  defined  as  the  plane  formed  by  the  vertical  axis  (as  described  by  grav¬ 
ity  direction)  and  the  vehicle  velocity  vector.  The  vehicle’s  forward  velocity,  V, 
is  assumed  constant,  and  the  model  tracks  the  vertical  position,  z,  the  transverse 
speed,  w,  along  with  the  pitch  angle,  6,  and  the  pitch  rate  q.  The  transverse  speed, 
w,  is  at  the  cavitator  and  is  defined  as  being  perpendicular  to  the  vehicle  axis. 
The  coordinate  system  is  attached  to  the  moving  vehicle  and  its  orientation  along 
with  the  speed  direction  definitions  are  shown  in  Figure  2.2.  The  control  action  is 
made  up  of  the  cavitator  and  fin  deflection  angles,  5C  and  5e,  respectively.  The  lift 
force  generated  by  the  control  elements  are  approximated  to  be  linearly  related  to 
their  angle  of  attack  with  respect  to  the  flow.  The  variable  n  is  the  ratio  of  lift 
effectiveness  between  the  fins  and  cavitator,  and  m  represents  the  density  ratio  of 
the  body  with  respect  to  the  surrounding  water.  Along  with  some  additional  small 
angle  assumptions,  the  equations  of  motion  are  represented  as  shown  in  Eqs.  (2.1) 
-  (2.4).  The  first  matrix  term  in  Eq.  (2.1)  includes  effects  due  to  body  orientation 
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and  the  coordinate  system.  The  second  term  incorporates  the  effect  of  the  control 
action.  The  third  matrix  includes  the  effect  of  gravity.  The  system  is  linear  with  the 
exception  of  the  planing  force  term,  Fp,  and  the  fourth  matrix  takes  into  account 
this  nonlinear  and  non-smooth  planing  force  that  occurs  when  the  body  contacts 
the  cavity  wall.  In  this  representation,  the  planing  force  is  taken  as  a  point  force 
applied  at  the  end  of  the  vehicle. 
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The  body  is  represented  as  a  conical  section  followed  by  a  cylindrical  section, 


as  shown  in  Figure  1.1.  In  this  formulation,  the  cavity  is  approximated  as  a  cylinder 
with  a  radius  equal  to  the  cavity  radius  at  the  rear  of  the  vehicle,  and  the  planing 
is  approximated  as  a  cylindrical  body  planing  on  a  cylindrical  cavity  surface.  The 


cavity  radius  at  the  rear  of  the  vehicle  is  calculated  by  using  a  formulation  that 
depends  only  on  the  cavitation  number,  a ,  and  the  cavitator  radius.  The  cavitation 
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number  is  a  non-dimensional  term  that  determines  the  extent  of  the  cavitation,  and 
this  number  is  defined  according  to  Eq.  (2.5).  In  this  expression,  p  represents  the 
fluid  density,  V  is  the  vehicle  velocity,  p^  is  the  ambient  fluid  pressure,  and  pc  is  the 
cavity  pressure.  The  model  is  run  with  a  vaporous  cavity,  and  since  the  fluid  and 
cavitator  parameters  are  constant,  the  cavity  radius  which  is  considered  to  depend 
only  on  the  forward  speed  of  the  vehicle  which  (within  each  simulation  run),  is 
constant.  In  this  form,  the  influence  of  the  cavitator  angle  of  attack  on  the  cavity 
radius  are  neglected. 


Poo  ~  Pc 

0.5  pV2 


(2.5) 


The  planing  force  function  is  shown  in  Eqs.  (2.3)-(2.4).  The  planing  force  is  a 


non-smooth  function  since  there  is  no  force  when  the  vehicle  is  not  in  contact  with 


the  cavity  wall.  The  quantity  h' ,  which  represents  the  immersion  depth,  is  the  source 
of  the  non-smoothness.  The  Rc  is  the  cavity  closure  rate  (at  the  rear  of  the  vehicle), 
which  is  included  as  a  correction  factor  to  the  angle  of  immersion,  a.  By  using 
the  cylinder-cylinder  assumption,  the  planing  force  can  be  determined  as  a  function 
of  the  angle  of  the  flow  velocity  with  respect  to  the  body.  Since  forward  speed  is 
constant,  this  angle  can  be  determined  from  the  transverse  speed.  A  representation 
of  the  planing  force  with  respect  to  the  transverse  speed  is  shown  in  Figure  2.3. 

This  system  can  show  unstable  behavior  when  no  control  is  present  [8],  [27]. 
By  using  a  simple  linear  feedback  controller,  stable  limit  cycles  can  be  observed  [8] , 
[27].  The  control  law  is  shown  in  Eqs.  (2.6),  and  only  the  cavitator  deflection  angle 
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Figure  2.3:  Original  planing  force  versus  transverse  speed. 


is  varied  while  holding  the  fin  deflection  constant  at  0.0  radians.  The  transverse 
velocity,  w  is  not  utilized  in  the  feedback  formulation  since  it  is  difficult  to  measure 
in  an  actual  application.  A  simulation  was  run  by  using  the  following  parameters: 
g  =  9.81  m/s2,  m  —  2,  Rn  =  0.0191  m,  R  —  0.0508  m,  L  —  1.8  m,  V  —  70.9740 
m/s,  a  =  0.0335,  n  =  0.5,  and  Cx o  =  0.82.  The  vehicle  parameters  are  chosen  to 
match  those  used  in  the  previous  literature  [8,  29,  30].  The  results  obtained  are 
presented  in  Figure  2.4.  Since  the  model  assumes  constant  forward  speed  (with 
respect  to  the  vehicle  axis),  the  transverse  speed,  w,  defines  the  velocity  angle  with 
respect  to  the  body.  Limit-cycle  motions  about  the  w  =  0  axis  indicate  the  potential 
for  two-sided  tailslap  behavior  where  the  vehicle  contacts  both  the  top  and  bottom 
surfaces  during  the  oscillatory  motion.  The  two-sided  planing  can  be  confirmed  for 
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Figure  2.4:  Simulation  run  demonstrating  limit-cycle  behavior  with  original  cavity 
model  and  planing  force  formulation;  controlled  case. 


this  simulation  run  by  tracing  the  vehicle  and  cavity  motions. 


Se  =  0 

Sc  =  15z  —  300  —  ,3q  (2.6) 


2.2  Integration  of  Numeric  Cavity  Model 

The  first  item  modified  is  the  cavity  model.  The  original  model  uses  a  closed- 
form  solution  for  the  cavity  radius.  The  numeric  cavity  model  described  in  references 
[19,  42]  is  then  implemented.  The  cavity  model  utilizes  an  iterative  potential  flow 
solver.  The  cavitation  number  and  cavitator  shape  (disc)  is  input  into  the  numeric 
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model  and  the  entire  non-dimensional  (normalized  with  respect  to  cavitator  diam¬ 
eter)  cavity  shape  for  axis-symmetric  flow  (no  angle  of  attack  for  the  cavitator)  is 
predicted. 

The  cavitator  angle  of  attack  also  affects  the  cavity  shape,  and  since  the  cavity 
model  only  provides  the  shape  of  the  cavity  for  an  axi-symmetric  flow,  the  angle  of 
attack  effects  are  included  as  a  refinement  to  the  axi-symmetric  data.  This  refined 
cavity  shape  is  calculated  as  a  shift  of  the  axi-symmetric  data  produced  by  the  cavity 
model.  The  shift  factor  is  a  term  derived  from  Logvinovich  [31],  and  this  factor  is 
based  on  the  principle  that  the  momentum  created  by  the  cavitator  lift  applies  an 
equal  and  opposite  momentum  upon  the  wake.  The  shift  can  be  expressed  as  in 
Eq.  (2.7),  where  Wy  expresses  the  lift  force  along  the  transverse  direction  (with 
respect  to  the  flow),  and  R  represents  unshifted  the  cavity  radius.  The  transverse 
lift  force  Wy  can  then  be  expressed  as  in  Eq.  (2.8).  In  this  expression,  a  represents 
the  cavitator  angle  of  attack  (with  respect  to  the  flow),  Dn  represents  the  cavitator 
diameter,  and  Cd  represents  the  coefficient  of  drag  for  axi-symmetric  flow. 


shift{x )  = 


W„ 


rx  ds 


npV2  Jo  R(s )2 
1  D2 

Wy  =  sin(a)  ■  cos(a)  -pn  —^C d 


(2.7) 

(2.8) 


When  normalized  with  respect  to  Dn,  the  shift  can  be  expressed  as  in  Eq.  (2.9). 
Here  the  '  designation  refers  to  the  non-dimensional  form,  so  that  shift'  =  shift/ Dn, 
R'  =  R/Dn,  and  x'  =  x/Dn.  The  cavity  model  generates  the  non-dimensional 
axi-symmetric  cavity  profile,  R'{x'),  as  well  as  the  coefficient  of  drag,  Cd-  The 
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Figure  2.5:  Non-dimensional  cavity  shape  with  angle  of  attack  effects. 

non-dimensional  cavity  profile  for  non-axi-symmetric  flow  can  then  be  expressed  by 
adding  the  shift,  as  in  Eq.  (2.10).  The  inclusion  of  the  cavity  shift  introduces  non¬ 
symmetry  and  cavity  radius  dependence  on  the  states  and  control.  An  example  of 
the  cavity  shift  is  shown  in  Figure  2.5.  The  flow  direction  is  taken  along  the  positive 
x-axis  and  an  angle  of  attack  of  10  degrees  is  assumed. 

Cn  rx'  ds 

shift1  (x )  =  — — sin(a)cos(a )  /  (2.9) 

8  Jo  K (s)z 

R'shifted(x )  =  R\x)  +  shift' (x)  (2.10) 

For  the  following  simulation  run,  the  original  planing  force  formulation  is  re¬ 
tained.  The  original  model  assumed  a  cylindrical  cavity  shape;  so  as  before,  only 
the  cavity  radius  at  the  rear  of  the  vehicle  (at  distance  of  L)  is  considered.  However, 
the  entire  cavity  profile  from  the  numeric  model  (as  well  as  the  Cd)  is  needed  to 
calculate  the  shift  associated  with  the  angle  of  attack.  The  cavitator  angle  of  attack 
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is  dependent  on  both  the  states  (which  determines  the  angle  of  attack  of  the  body 
with  respect  to  the  flow)  and  the  control  (cavitator  angle  of  attack  with  respect 
to  the  body).  The  cavity  radius  is  no  longer  constant  and  varies  throughout  the 
simulation.  Additionally,  since  the  shift  is  non-symmetric,  the  cavity  shape  is  not 
symmetric,  and  a  distinction  must  be  made  between  planing  on  the  top  surface  or 
bottom  surface. 

2.2.1  Simulation  Runs  with  Numeric  Cavity  Model 

The  simulations  were  run  with  the  same  parameters  as  in  the  initial  simulation; 
that  is,  g  =  9.81  m/s2,  m  —  2  kg,  Rn  =  0.0191  m,  R  =  0.0508  m,  L  =  1.8  m, 
V  =  70.9740  m/s,  a  =  0.0335,  n  =  0.5,  and  Cx 0  =  0.82.  The  results  are  presented  in 
Figure  2.6  and  a  qualitative  difference  in  the  system  behavior  is  evident.  Limit  cycles 
are  not  present.  This  is  partially  related  to  the  difference  in  the  nominal  (unshifted) 
cavity  radius  calculated  from  the  numeric  model.  A  comparison  of  the  cavity  radius 
versus  cavitation  number  is  shown  in  Figure  2.7.  The  numeric  model  provides  a 
larger  cavity  radii  than  the  original  model.  The  larger  radius  creates  a  larger  cavity 
to  body  clearance  at  a  given  speed  (cavitation  number).  The  larger  cavity  to  wall 
clearance  is  enough  to  transition  the  system  out  of  the  limit  cycle  region  when  the 
original  run  parameters  are  used.  It  should  be  noted  that  the  new  cavity  model 
does  not  eliminate  limit-cycle  behavior  in  all  scenarios.  If  the  speed  were  decreased, 
or  the  body  radius  increased  (to  cause  lower  body  to  wall  clearances),  limit  cycles 
can  be  observed.  Limit-cycle  observations  will  be  discussed  in  a  subsequent  section. 
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Time  (s)  Time  (s) 


Figure  2.6:  Simulation  run  with  numeric  cavity  model  and  original  planing  force 
formulation  with  linear  feedback. 

However,  the  qualitative  difference  in  behavior  when  using  the  same  run  parameters, 
does  show  the  significance  of  the  cavity  model.  Anglc-of-attack  effects  can  be  seen  in 
Figure  2.8,  wherein  the  effective  cavity  radius  (depending  on  top  or  bottom  planing) 
observed  during  the  simulation  is  shown. 


2.3  Planing  Force  Model 

The  second  item  to  be  refined  is  the  planing  force  formulation.  Three  plan¬ 
ing  force  models  were  investigated.  The  models  were  developed  by  Hassan  [14], 
Logvinovich  [32],  and  Paryshev  [37].  Dzielski  presents  an  interpretation  of  the  three 
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Rc  vs.  captation  number 


Figure  2.7:  Cavity  radius  versus  cavitation  number  as  obtained  from  the  cavity 
model  used  in  earlier  studies  [8,  29,  30]  and  the  numerical  model.  Planing  force 
direction  convention  taken  with  respect  to  the  positive  transverse  velocity  direction 

w. 
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Cavity  Radius  vs.  Time 


Figure  2.8:  Effective  cavity  radius  (depending  on  top/bottom  planing)  versus  time 


for  simulation  run  by  using  numeric  cavity  model  with  angle  of  attack  effect. 
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approaches  in  relation  to  a  cylinder  planing  on  a  cylindrical  surface  [7].  Given  the 
definition  of  terms  provided  in  Figure  2.9,  the  apparent  mass  can  be  calculated  ac¬ 
cording  to  Eq.  (2.11),  where  A  —  R  —  r.  The  Hassan  and  Logvinovich  models  make 
use  of  m*L  as  the  apparent  mass,  while  the  Paryshev  model  makes  use  of  mP.  The 
derivative  with  respect  to  h  is  then  taken,  and  Eq.  (2.12)  is  used  in  integrating 
the  planing  force.  For  a  cylinder  planing  on  a  cylinder,  Eq.  (2.13)  represents  the 
planing  force.  A  closed-form  solution  for  the  integral  can  be  found  for  the  Hassan 
and  Paryshev  models.  The  Hassan  model  simplifies  to  the  planing  force  model  used 
in  the  previous  works  (as  presented  in  the  previous  section).  The  Paryshev  solution 
for  the  cylinder-cylinder  case  is  given  in  Eq.  (2.14),  where  XP  is  the  force  centroid. 
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Figure  2.9:  Diagram  of  a  cylinder  planing  on  a  cylindrical  surface. 

The  three  methods  provide  different  planing  forces.  A  comparison  of  the  plan¬ 
ing  force  results  (in  relation  to  w  in  the  system  model)  is  shown  in  Figure  2.10.  The 
Logvinovich  results  were  calculated  by  using  a  numerical  integration.  Both  Pary- 
shev  and  Logvinovich  produce  higher  planing  forces  than  the  planing  formulation 
used  in  the  original  model.  The  Paryshev  method  was  chosen  for  integration  into 
the  model  since  it  was  shown  to  better  fit  experimental  results  for  both  static  and 
dynamic  planing  [7]. 

All  of  these  planing  force  representations  are  for  steady  planing  (no  body 
acceleration  or  velocity  into  or  out  of  the  fluid).  A  derivation  of  the  Paryshev 
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Planing  force  vs  vertical  speed 


Figure  2.10:  Comparison  of  planing  force  models  for  cylindrical  cavity  shapes. 

planing  force  model  as  applied  in  this  dissertation  is  presented  next. 

2.3.1  Paryshev  Planing  Force 

Paryshev  defines  the  forces  on  an  expanding  cylinder  planing  on  a  cylindrical 
cavity.  The  force  per  unit  length  is  given  by  the  rate  of  change  of  momentum  of 
the  fluid  displaced  by  the  planing  cylinder.  This  is  shown  in  Eq.  (2.15),  where 
m*  represents  the  apparent  mass  due  to  the  planing,  and  mR  is  the  apparent  mass 
due  to  the  expansion  of  the  cylinder.  For  the  supercavitating  vehicle  dynamics,  the 
contribution  of  the  m*R  can  be  ignored  since  the  radius  of  the  body  does  not  change. 
The  expression  for  the  force  can  then  be  expanded  as  shown  in  Eq.  (2.16).  If  only 
steady  planing  is  considered  (there  is  no  acceleration  into  or  out  of  the  fluid)  the 
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second  term  involving  ^  can  be  dropped. 
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The  expression  for  the  apparent  mass  m*  is  given  in  Eq.  (2.17),  and  this 
expression  is  a  function  of  the  immersion  depth  h.  The  gap  A  is  defined  as  A  =  R—r. 
The  rate  of  change  of  the  apparent  mass  m*  can  then  be  described  as  shown  in  Eq. 
(2.18).  Again,  since  the  radius  of  the  body  does  not  change,  the  gap  is  constant 
with  time,  so  ^  =  0.  The  term,  -jjr,  can  be  expressed  as  in  Eq.  (2.19).  The  total 
planing  force  can  then  be  determined  by  integrating  over  the  entire  wetted  area  as 
represented  as  in  Eq.  (2.20),  with  Vy  =  h. 
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For  the  case  of  a  cylindrical  cavity,  the  planing  force  can  be  integrated  over  the 
planing  area  and  the  obtained  solution  is  given  in  Eq.  (2.21).  In  these  equations, 
ho  is  the  immersion  depth  at  the  aft  of  the  vehicle,  as  shown  in  Figure  2.9.  The 
immersion  rate  term  h  can  be  expressed  as  Vsin(a).  The  force  centroid  can  also  be 
calculated,  as  shown  in  Eq.  (2.22),  with  Xp  measured  from  the  aft  of  the  vehicle. 
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Again,  the  Paryshev  planing  force  representation  was  chosen,  as  it  has  been  shown 


to  provide  a  better  fit  to  experimental  planing  force  data  [7]. 
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2.4  Simulation  Results  with  Numeric  Cavity  Model  and  Paryshev 
Planing  Force  Formulation 

By  utilizing  the  new  cavity  and  planing  force  representations  additional  simu¬ 
lations  can  be  carried  out.  The  same  parameters  as  used  in  the  previous  simulation 
runs  are  utilized,  g  =  9.81  m/s2,  rn  —  2  kg,  Rn  =  0.0191  m,  R  =  0.0508  m,  L  =  1.8 
m,  V  =  70.9740  m/s,  a  =  0.0335,  n  =  0.5,  and  Cx0  =  0.82.  The  same  linear  feed¬ 
back  control  law  is  implemented  as  given  by  Eqs.  (2.6),  where  only  the  cavitator  is 
actuated  and  the  fins  are  assumed  to  be  passive.  The  same  parameters  are  utilized 
to  present  a  direct  comparison  with  the  previous  simulation  run  shown  in  Figure 
2.4,  where  two-sided  tailslap  behavior  is  observed. 

The  simulation  results  are  presented  in  Figure  2.11.  At  steady  state,  limit- 
cycle  oscillations  are  observed.  The  steady-state  oscillations  in  w  do  not  cross  the 
w  =  0  axis.  Assuming  a  small  cavity  shift,  this  indicates  that  there  is  no  two- 
sided  tailslap  behavior  (oscillations  where  the  cavity  boundaries  are  crossed  on  both 
sides).  The  cavity  and  supercavitating  body  motions  can  be  tracked  to  confirm  that 
the  oscillations  are  indeed  only  about  the  bottom  planing  surface.  The  lack  of  two- 
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Figure  2.11:  Results  obtained  with  cylindrical  planing  force  formulation  and  linear 
feedback. 

sided  tailslapping  is  due  to  the  different  cavity  and  planing  force  models  used  here. 
As  discussed  in  the  previous  section,  the  cavity  model  used  in  this  effort  provides 
a  slightly  larger  nominal  (un-shifted)  cavity  radius  than  that  used  in  the  previous 
studies  [8,  29,  30]  (Figure  2.7),  and  the  increased  cavity  to  body  clearance  is  partially 
responsible  for  inhibiting  the  two-sided  tailslap  behavior. 
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Figure  2.12:  Results  obtained  with  cylindrical  planing  force  formulation  and  double 
linear  feedback  in  the  presence  of  downwash  effects. 

2.4.1  Cavity  Shift  Effect 

The  difference  in  system  behavior  is  not  only  due  to  the  change  in  nominal 
cavity  radius.  With  this  cavity  formulation,  the  cavity  shape  is  not  constant  as  in 
previous  studies  [8,  29,  30].  The  effect  of  the  cavity  shift  can  be  noted  with  a  stabi¬ 
lized  system.  By  using  the  previous  run  parameters,  the  system  can  be  stabilized, 
as  shown  in  Figure  2.12,  by  doubling  the  values  of  the  control  law  coefficients  given 
in  Eqs.  2.6.  However,  even  with  this  control,  if  the  cavity  shift  effect  is  removed, 
the  system  can  again  be  seen  to  exhibit  oscillations. 
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Planing  force  versus  vertical  speed 


Figure  2.13:  Normalized  planing  force  versus  vertical  speed  w,  with  cavitator  actu¬ 
ation  neglected,  for  cylindrical  planing  force  model. 

The  cavity  shift  is  dependent  on  the  angle  of  attack  of  the  cavitator  with 
respect  to  the  flow.  The  cavitator  angle  of  attack  is  in  part  due  to  the  vehicle  orien¬ 
tation  (with  respect  to  the  flow).  As  the  vehicle  beings  to  plane,  the  contribution  of 
the  vehicle  orientation  on  the  cavitator  angle  of  attack  will  tend  to  shift  the  cavity 
away  from  the  body.  A  representation  of  the  normalized  planing  force  versus  verti¬ 
cal  speed  w  is  shown  in  Figure  2.13.  Here,  the  planing  forces  with  and  without  the 
shift  effect  are  shown.  The  cavitator  actuation  angle  is  fixed  at  Sc  =  0  to  remove 
contributions  due  to  the  actuation  angle.  As  shown,  the  cavity  shift  due  to  body 
orientation,  can  be  considered  as  creating  a  slightly  larger  cavity  radius,  which  in 
some  cases,  may  delay  the  onset  of  oscillations. 
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Chapter  3 


Vehicle  Dynamics  with  Shifted  Non-Cylindrical  Cavities 

Since  the  cavity  model  produces  an  entire  cavity  shape,  the  planing  force  can 
be  calculated  on  the  basis  of  the  entire  profile  rather  than  on  the  basis  of  a  cylin¬ 
drical  shape  assumption,  which  is  fairly  inaccurate  in  determining  the  planing  area, 
particularly  during  high  immersion.  With  a  cylindrical  cavity  shape  assumption, 
the  wetted  area  is  generally  over-predicted  as  shown  in  Figure  3.1.  The  horizontal 
line  represents  the  assumed  cavity  shape,  with  the  cylindrical  cavity  assumption. 
The  curved  line  represents  the  actual  cavity  shape.  The  error  in  cavity  radius  made 
by  using  the  cylindrical  assumption  during  vehicle  system  simulations  can  be  quite 
high  [35]. 

3.1  Non-cylindrical  Cavity  Planing 

To  incorporate  the  non-cylindrical  profile,  as  shown  in  Figure  3.2,  the  cavity 
can  be  treated  as  made  up  of  several  short  cylindrical  sections.  During  each  time 
step,  the  entire  shifted  cavity  profile  is  determined.  Areas  of  interference  with  the 
vehicle  body  are  then  determined.  No  planing  is  assumed  to  occur  along  the  conical 
forebody.  The  planing  area  is  then  re-meshed  by  using  linear  interpolation  between 
the  cavity  points  to  refine  the  region  of  planing. 

The  cavity  contraction  rate  must  also  be  considered  for  the  planing  force. 
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Figure  3.1:  Illustration  for  cylindrical  cavity  assumption. 

In  the  planing  force  formulation,  only  the  angle  of  planing  a  is  used  to  generate 
the  immersion  rate  term  h.  So  for  planing  forces,  the  contraction  rate  can  be 
expressed  as  an  augmentation  to  a  (similar  to  the  role  of  Rc  in  reference  [8]).  If 
independent  expansion  is  assumed,  as  in  the  cavity  model  presented  in  reference 
[31],  the  contraction  rate  can  be  determined  by  the  cavity  slope  and  velocity.  In 
the  formulation  for  non-cylindrical  planing,  the  cavity  slope  and  the  body  slope 
along  each  integration  section  are  used  to  calculate  an  effective  a  (see  Eq.  (3.1)), 
which  generates  an  appropriate  h  (for  each  section)  that  accounts  for  both  the  body 
immersion  rate  and  the  cavity  contraction.  The  planing  force  is  then  numerically 
integrated  across  all  sections  of  interference. 

aCorrected{x)  =  a  +  tan-1  (SR(x)/S(x))  (3.1) 

In  the  cylindrical  planing  formulation,  since  the  cavity  shape  is  assumed,  the 
cavity  shift  only  affects  when  the  vehicle  begins  to  plane  and  it  does  not  affect 
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Figure  3.2:  Cavity  approximated  by  a  series  of  short  cylindrical  sections. 

the  planing  force  function  once  planing  has  been  initiated.  By  contrast,  in  the 
non-cylindrical  planing  force  formulation,  the  shift  not  only  affects  when  the  vehicle 
begins  to  plane,  but  it  also  changes  the  profile  of  the  cavity  and  therefore  the  planing 
force  function  once  planing  has  been  initiated.  A  depiction  of  the  planing  force 
dependence  is  shown  in  Figure  3.3.  The  control  as  well  as  the  states  (which  determine 
body  orientation),  determine  the  cavitator  angle  of  attack.  This  in  turn  determines 
the  shift  and  the  cavity  profile.  The  profile  as  well  as  the  body  orientation  are  used 
in  determining  whether  the  body  is  in  contact  with  the  cavity  or  not.  If  there  is 
contact,  the  profile  and  the  body  orientation  are  utilized  to  integrate  the  planing 
force.  In  previous  research  [29],  wherein  cavity  shift  effects  were  not  considered, 
the  vehicle  dynamics  could  be  modeled  by  using  a  switching  system,  with  a  set  of 
equations  governing  the  system  when  planing  occurs  and  another  set  of  equations 
when  planing  does  not  occur.  In  these  cases,  there  are  defined  state  dependent 
switching  boundaries  and  defined  switched  dynamics.  However,  when  one  considers 
the  cavity  shift  as  well  as  the  non-cylindrical  planing,  the  system  has  both  a  state 
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Figure  3.3:  Planing  force  variation  for  the  non-cylindrical  planing  force  model. 


and  control  dependent  switching  boundary,  as  well  as  a  state  and  control  dependent 
switched  dynamics. 

The  planing  force  function  is  represented  in  Figure  3.4.  Here,  the  planing  force 
variation  is  shown  with  respect  to  the  vertical  speed  of  the  body  w,  as  well  as  the 
cavitator  actuation  angle  Sc.  Cavity  shift  effects  are  incorporated,  the  cavity  shift 
being  a  function  of  both  the  body  orientation  (which  for  a  constant  forward  speed 
becomes  solely  a  function  of  w),  and  the  cavitator  actuation  angle.  The  onset  of 
planing  in  relation  to  w  changes  with  respect  to  the  actuation  angle.  Once  planing 
is  initiated,  the  planing  force  can  also  be  seen  to  vary  with  respect  to  the  actuation 
angle. 
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Figure  3.4:  Normalized  planing  force  versus  vertical  speed  w,  and  cavitator  actuation 
angle  Sc,  for  non-cylindrical  planing  force  model. 

3.2  Non-cylindrical  Cavity  Simulation  Results 

Results,  which  were  obtained  by  using  the  same  vehicle  parameters  as  before 
and  the  feedback  control  law  given  by  Eqs.  (2.6),  are  presented  in  Figure  3.5. 
Oscillations  about  a  single  planing  surface  are  observed.  In  prior  research  with 
cylindrical  cavity  models  [29],  it  was  found  that  when  the  cavitation  number  is 
lowered  from  a  =  0.0335  to  a  =  0.025  (increased  speed  and  larger  cavities),  the 
system  can  demonstrate  stable  equilibrium  solutions.  A  similar  remark  can  be 
made  about  the  non-cylindrical  planing  force  case.  Results  obtained  for  cr  =  0.025 
are  shown  in  Figure  3.6.  Here,  the  system  appears  to  stabilize  to  an  equilibrium 
solution.  Additionally,  when  the  cavitation  number  is  increased  from  a  =  0.0335 
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Figure  3.5:  System  response  results  obtained  with  non-cylindrical  planing  force 
formulation  and  linear  feedback.  Cavitation  number  is  cr  =  0.0335. 

to  a  =  0.037  (decreased  speed  and  smaller  cavities),  the  vehicle  is  seen  to  exhibit 
two-sided  tailslap  motions,  as  shown  in  Figure  3.7. 

In  prior  efforts  in  the  group,  a  Hopf  bifurcation  of  an  equilibrium  solution  and 
a  period-doubling  bifurcation  were  identified  [29].  In  order  to  determine  whether 
these  bifurcations  also  occur  in  the  system  with  non-cylindrical  and  non-symmetrical 
cavities,  an  attempt  was  made  to  numerically  determine  the  equilibrium  solutions. 
However,  even  for  cases,  where  the  time  domain  results  suggest  the  presence  of  a 
stable  equilibrium  solution,  these  equilibrium  solutions  could  not  be  found.  This 
is  suspected  to  be  due  to  the  discretization  of  the  sections  used  for  the  numerical 
integration  along  the  cavity  profile. 


42 


Figure  3.6:  System  response  results  obtained  with  non-cylindrical  planing  force 
formulation  and  linear  feedback.  Cavitation  number  is  a  =  0.025. 

3.3  Smoothened  Non-Cylindrical  Cavities  and  Equilibrium  Points 

Even  with  the  re-meshed  area  of  planing,  as  outlined  in  the  previous  section, 
only  a  discrete  estimate  of  the  planing  area  (and  subsequent  planing  force)  is  ob¬ 
tained.  To  improve  this  computation,  a  smoothened  version  of  the  cavity  can  be 
formed  by  carrying  out  a  cubic  spline  interpolation  of  the  cavity  coordinates  gener¬ 
ated  by  the  numerical  cavity  model.  This  creates  a  third-order  piecewise  polynomial 
representation  of  the  cavity  profile.  With  a  continuous  expression  for  the  entire  cav¬ 
ity,  when  the  body  is  immersed,  the  planing  area  and  planing  force  can  be  expressed 
as  smooth  functions.  It  should  be  noted  that  unlike  the  smoothing  used  in  reference 
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Figure  3.7:  System  response  results  obtained  with  non-cylindrical  planing  force 
formulation  and  linear  feedback.  Cavitation  number  is  a  =  0.037. 

[28],  here,  the  transition  to  planing  is  still  non-smooth.  Planing  forces  determined 
with  the  smoothened  cavity  model  and  the  discrete  section  model  are  compared  in 
Figure  3.8,  for  zero  cavitator  actuation  angle;  that  is,  Sc  =  0.  There  are  no  dis¬ 
cernible  differences,  with  the  remark  that  the  chattering  type  behavior  is  eliminated 
with  the  splincd  cavity  approximation.  Similar  results  were  also  obtained  for  other 
cavitator  actuation  angles.  The  time  domain  results  obtained  with  the  splined  cavity 
are  also  matched  with  those  obtained  earlier  for  the  discrete  cavity  sections. 

The  main  advantage  of  using  a  splined  cavity  profile  is  that  equilibrium  solu¬ 
tions  can  be  obtained.  A  branch  of  equilibrium  points  related  to  stable  operations  is 
shown  in  Figure  3.9.  In  this  figure,  the  stars  represent  stable  equilibrium  points  and 
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Planing  Forces  Using  Linear  Interpolation  and  Splined  Cavity,  6C  =  0 


Figure  3.8:  Normalized  planing  force  versus  vertical  speed,  w,  for  both  the  re-meshed 
discrete  sections  and  splined  cavity  shapes,  when  5C  =  0. 

the  circles  represent  unstable  equilibrium  points.  It  can  be  seen  that  at  a  cavitation 
number  between  0.0313  and  0.0315,  there  is  a  transition  from  a  stable  state  to  an 
unstable  state  on  this  branch  of  equilibrium  points.  Here  the  equilibrium  solutions 
are  plotted  as  the  l2  norms  of  the  corresponding  states.  Thus,  each  solution  can  be 
represented  as  a  scalar  value,  which  corresponds  to  the  magnitude  of  an  equilibrium 
solution;  this  is  done  for  easier  visualization. 

Through  the  time  domain  results,  stable  equilibrium  and  stable  limit-cycle 
behavior  are  identified.  In  Figure  3.10,  a  graph  of  the  projection  of  the  steady-state 
behavior  of  the  system  is  shown  on  the  w  -  q  plane  with  respect  to  the  cavitation 
number.  As  expected,  the  transition  from  stable  to  unstable  equilibrium  points 
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Figure  3.9:  l 2  norm  of  the  equilibrium  points  versus  cavitation  number. 

corresponds  to  the  onset  of  limit-cycle  behavior  and  therefore  suggesting  a  Hopf  bi¬ 
furcation.  By  numerically  evaluating  the  eigenvalues  of  the  linearized  system  about 
the  equilibrium  points,  the  system  is  seen  to  satisfy  the  transversality  condition  re¬ 
quired  for  a  Hopf  bifurcation.  This  finding  parallels  the  Hopf  bifurcation  found  in 
the  prior  work  with  a  cylindrical  cavity  formulation  [28,  29].  However,  the  instability 
occurs  at  a  different  cavitation  number  in  the  present  case.  As  the  cavitation  num¬ 
ber  is  further  increased,  an  abrupt  change  in  the  character  of  the  steady-state  limit 
cycles  is  observed,  as  shown  in  Figure  3.11.  At  a  cavitation  number  between  0.036 
and  0.03625,  the  steady-state  limit-cycle  motion  grows  dramatically.  From  the  time 
domain  responses  of  the  vehicle,  it  is  found  that  this  abrupt  change  corresponds  to 
the  case  when  the  limit  cycles  first  begin  to  plane  about  both  surfaces  at  steady 
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Figure  3.10:  Projection  of  steady-state  behavior  of  system  in  the  w-q  plane  versus 
cavitation  number  a. 

state.  It  is  believed  that  this  additional  impact  with  the  second  cavity  surface  is 
related  to  the  post-Hopf  bifurcation  behavior  discussed  in  reference  [29]. 

3.3.1  System  Dynamics  With  Washout  Filter 

In  the  previous  work  reported  in  reference  [29],  a  washout  filter  was  used  to 
delay  the  onset  of  the  Hopf  bifurcation.  This  approach  is  desirable  since  it  preserves 
the  equilibrium  solutions  of  the  original  system.  The  washout  filter  generates  an 
additional  state  variable  p  whose  dynamics  are  described  by  Eq.  3.2,  where  d  >  0 
is  a  constant  between  0  and  1,  kp  is  the  dynamic  feedback  coefficient,  and  q  is  the 
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Figure  3.11:  Projection  of  steady-state  behavior  of  system  in  the  w-q  plane  versus 
cavitation  number  a,  showing  two-sided  tailslap  behavior. 

pitch  rate.  The  original  feedback  law  given  by  Eqs.  (2.6)  is  then  augmented  by 
adding  <5c2,  as  given  in  Eq.  (3.3). 


p  =  q  —  dp 

(3.2) 

kp(q  -  dp) 

(3.3) 

For  an  equilibrium  solution  q  —  dp  =  0  and  the  system  simplifies  back  to  the 
original  system  without  the  filter.  Equilibrium  solutions  obtained  for  d  =  0.5  and 
kp  =  —3,  are  shown  in  Figure  3.12.  These  equilibrium  solutions  are  the  same  as  those 
obtained  in  the  original  system  and  shown  in  Figure  3.9,  but  here,  these  equilibrium 
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Figure  3.12:  l 2  norm  of  the  equilibrium  points  versus  cavitation  number  with 
washout  filter. 

points  remain  stable  and  there  are  no  limit-cycle  motions  within  the  considered 
range  of  the  cavitation  number.  A  time  domain  response  of  the  system  with  the 
washout  filter  obtained  at  a  =  0.034  is  shown  in  Figure  3.13.  The  steady  state 
here  corresponds  to  a  stable  equilibrium  position,  while  in  the  original  system,  this 
cavitation  number  is  associated  with  limit-cycle  motions  in  the  post-Hopf  bifurcation 
regime. 

3.3.2  Fin  Input  Based  Stabilization  Inside  Cavity 

Stabilization  inside  the  cavity  may  also  be  beneficial  for  straight  line  flight. 
Since  there  are  no  planing  forces  while  the  vehicle  is  completely  inside  the  cavity, 
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Figure  3.13:  Results  obtained  with  non-cylindrical  planing  force  formulation  and 
washout  filter.  Cavitation  number  is  a  =  0.034. 

the  fins  need  to  be  utilized  to  stabilize  the  rear  of  the  vehicle.  Through  the  previous 
simulations,  it  can  be  seen  that  a  fin  deflection  of  Se  =  0  is  not  enough  to  support 
the  rear  of  vehicle.  Instead,  a  passive  hn  deflection  of  Se  =  0.1  is  utilized  along 
with  the  cavitator  control  input  described  by  Eqs.  2.6.  The  hn  deflection  angle 
was  chosen  by  determining  the  hn  force,  under  static  conditions,  that  is  required 
to  support  the  rear  of  the  vehicle.  The  vehicle  system  is  found  to  have  stable 
equilibrium  solutions  inside  the  cavity.  The  equilibrium  points  are  plotted  in  Figure 
3.14.  Stable  equilibrium  points  exist  past  the  region  where  Hopf  bifurcation  occurs 
in  the  original  system. 
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Figure  3.14:  l 2  norm  of  the  equilibrium  points  versus  cavitation  number  by  using 
linear  feedback  and  constant  fin  input. 

3.4  Logvinovich  Cavity  and  Related  Dynamics 

A  similar  analysis  was  also  performed  by  using  a  second  cavity  model.  A 
semi-empirical  closed-form  solution  for  cavity  shape  is  presented  in  reference  [31]. 
A  diagram  of  the  cavity  is  shown  in  Figure  3.15.  For  this  representation,  the  cavity 
radius  is  expressed  as  in  Eq.  (3.4).  Here,  Lk  and  Rk  occur  at  the  maximum  radius 
of  the  cavity.  The  terms  X\  and  R\  refer  to  an  arbitrary  point  along  the  cavity, 
and  x  is  a  correction  factor.  The  cavity  contour  is  defined  as  passing  through  the 
point  where  R  =  R\  at  time  t  —  0.  An  additional  derived  relationship  between  Rk 
and  Rn  is  shown  in  Eq.  (3.5).  In  this  expression,  cx o  =  0.82  and  the  constant  k  is 
approximately  unity  (for  the  cavitations  numbers  of  interest);  both  of  these  values 
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Figure  3.15:  Diagram  for  the  Logvinovich  Cavity  Model 

can  be  determined  experimentally.  Additionally,  if  x\  is  chosen  as  x\  =  2 Rn,  R± 
can  be  expressed  as  in  Eq.  (3.6).  Furthermore  from  Figure  3.15,  tk  =  L^/V,  and  t 
can  be  expressed  in  terms  of  the  distance  a;  as  t  =  (x  —  x\ )/V.  The  term  Lk  can 
in  turn  also  be  approximated  by  using  an  experimentally  derived  relationship  as  in 
Eq.  (3.7). 


R  —  Ri 


\ 


1  - 


'  rT 
,  Rl 


i  -  — 

tk 


2/X 


Rk\  2  _  CX o(l  +  cr) 
Rn  J  ka 

R\  =  1.92i?n 

1.92 


7  R' 


-3 


a 


(3.4) 

(3.5) 

(3.6) 

(3.7) 


With  a  correction  factor  of  x  =  0.85  (chosen  to  match  experimental  data), 
the  expression  for  cavity  radius  can  be  written  as  Eq.  (A.  14),  where  the  d  terms 
represent  the  associated  diameters. 
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dmax  =  dc\]  0.82(1  +  a)/a 

lm  =  dc/2(1.92/a  -  3) 
k\  =  1.92(.82(1  +  <j)/<j)~^ 

*  dc  dc)  j Ifji 

Rc(x)  =  dmax/2^/l-(l-kl)\l-k2\^  (3.8) 

A  comparison  of  cavity  models  is  presented  in  Figure  3.16.  Here,  three  cavity 
models  are  presented,  the  Logvinovich  model  described  in  this  section,  the  numeric 
cavity  model,  and  the  Dziclski  cavity  representation  presented  in  reference  [8]  (only 
used  for  approximating  of  cavity  radius  at  the  rear  of  the  vehicle).  For  the  vehicle 
length  considered  in  this  work  (L  =  1.8m),  the  Logvinovich  and  numeric  cavity 
models  produce  similar  results  at  the  tail  end  of  the  vehicle. 

The  Logvinovich  cavity  formulation  was  implemented  into  the  non-cylindrical 
planing  models  by  sampling  points  along  the  cavity  length.  The  simulations  provided 
qualitatively  similar  results  as  those  obtained  with  the  numeric  cavity  model.  The 
system  exhibits  stable  equilibrium  behavior  at  low  cavitation  numbers  (high-speed) 
and  transitions  into  limit-cycle  motion  as  the  cavitation  number  increases  and  the 
cavities  tighten.  A  diagram  of  the  equilibrium  solutions  is  presented  in  Figure  3.17. 
The  system  is  also  similarly  stabilized  by  using  a  washout  filter  or  fin  input.  A 
full  write-up  on  the  simulation  results  using  the  Logvinovich  cavity  is  presented  in 
reference  [34]. 

The  main  advantage  of  the  Logvinovich  cavity  representation  is  that  it  is  a 
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Figure  3.16:  Cavity  model  comparison  for  a  =  0.0335. 

closed-form  solution,  and  cavity  shapes  can  be  easily  and  rapidly  calculated.  This 
aspect  is  not  as  important  for  constant  speed  simulations  (as  presented  in  the  pre¬ 
vious  sections),  but  becomes  important  for  variable  speed  maneuvering  which  will 
be  presented  in  a  subsequent  chapter. 
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Figure  3.17:  l 2  norm  of  the  equilibrium  points  versus  cavitation  number  by  using 
the  Logvinovich  cavity  formulation. 
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Chapter  4 


Extensions  to  Vehicle  Dynamics  Models  and  Related  Dynamics 

In  this  chapter,  some  additional  model  considerations  are  presented  for  the 
supercavitating  vehicle  system.  The  modifications  and  the  physical  basis  behind  the 
modifications  are  discussed,  along  with  an  examination  of  the  resulting  dynamics. 

4.1  Partial  Cavitation 

Being  able  to  model  dynamics  for  partially  cavitating  vehicles  is  an  important 
step  towards  being  capable  of  modeling  full  vehicle  missions.  Partial  cavitation  oc¬ 
curs  when  operating  at  speeds  where  the  cavitation  bubble  does  not  entirely  envelop 
the  vehicle,  as  illustrated  in  Figure  4.1.  When  operating  in  these  ranges,  a  portion 
of  the  vehicle  (from  the  rear  forward)  is  fully  wetted.  Along  the  wetted  portions,  the 
vehicle  experiences  buoyancy  forces,  and  added  mass  effects  (caused  by  accelerating 
the  body  in  the  fluid).  The  planing  no  longer  occurs  at  the  rear,  but  occurs  where 
the  cavity  is  formed. 

For  the  partial  cavitation  case,  a  more  general  forebody  is  considered,  as  shown 
in  Figure  4.2.  A  truncated  cone  with  an  initial  radius  of  R\  is  considered.  This 
modification  was  made  because  the  numeric  cavity  models  used  to  predict  body-in¬ 
flow  partial  cavities  were  made  with  these  truncated  cone  forebodies. 
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Figure  4.2:  Vehicle  and  cavity  length  parameters. 


4.1.1  Buoyancy  and  Added  Mass  Effects 

If  the  dimensions  are  as  shown  in  Figure  4.2,  the  buoyancy  force  can  be  de¬ 
termined  as  shown  in  Eq.  (4.1),  where  rc  =  Lc  +  r\.  The  resulting  buoyancy 
moments  can  be  expressed  as  given  in  Eqs.  (4.2)-(4.4). 


pgirR2(L  -  Lc )  if  Lc  >  Lx 

Fbuoy  =  <  (4.1) 

[  pgn  (R2{L  -  M)  +  \  {LX  -  Lc)  otherwise 
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(4.2) 
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Mbuoy2  otherwise 
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(4.3) 

(4.4) 


The  added  mass  terms  are  approximated  as  the  added  mass  of  a  2-D  cylinder 
moving  through  fluid.  A  similar  approach  has  been  presented  in  reference  [31].  The 
obtained  force  and  moments  due  to  the  added  mass  effects  are  shown  in  Eqs.  (4.5)- 
(4.6).  The  added  mass  terms  are  defined  according  to  Eqs.  (4.7)-(4.9),  and  the  force 
and  moment  conventions  remain  the  same  as  with  the  previous  systems  with  the 
reference  being  the  nose  of  the  vehicle. 


Fam  ~  X22U!  + 

(4.5) 

Mam  =  X26W  +  X^q 
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(4.8) 
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4.1.2  Cavity  Model 


For  partial  cavitation,  the  cavity  prediction  model  must  now  account  for  the 
vehicle  body  which  is  now  within  the  flow.  Here,  a  partial  cavitating  numeric  cavity 
model  provided  in  reference  [44]  is  utilized.  This  numeric  model  can  be  used  to  solve 
for  the  steady-state  cavities  that  close  along  a  body  for  an  axi-symmetric  flow,  as 
shown  in  Figure  4.3. 

The  cavity  model  is  an  iterative  potential  flow  solver  similar  to  the  supercav- 
itating  numeric  cavity  model.  Cavity  shift  effects  due  to  cavitator  angle  of  attack 
are  again  approximated  by  using  the  method  shown  in  the  previous  chapter.  A 
limitation  of  this  model  is  that  it  is  only  able  to  converge  for  small  cavities  (high 
cavitation  numbers).  So  only  cavities  that  close  on  the  front  portions  of  the  body 
can  be  considered  in  this  fashion.  As  shown  in  Figure  4.1,  planing  occurs  along 
the  length  of  the  cavity,  and  where  the  vehicle  is  wetted,  added  mass  and  buoyancy 
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Figure  4.3:  System  response  results  from  numeric  partial  cavity  model  presented  in 
reference  [44], 

forces  are  introduced  into  the  system  dynamics. 

4.1.3  Simulation  Results 

The  partially  cavitating  vehicle  system  was  run  by  using  the  same  simulation 
parameters  as  used  in  the  previous  studies,  with  g  =  9.81  m/s2,  m  —  2  kg ,  Rn  = 
0.0191  m,  R  —  0.0508  m,  L  =  1.8  m ,  n  =  0.5,  and  Cx o  =  0.82.  Working  near  the 
largest  cavity  limit  of  the  partial  cavity  model  (for  the  given  vehicle  and  cavitator 
parameters),  a  cavitation  number  of  a  =  0.066925  is  considered.  This  generates 
a  cavity  of  approximately  0.7640  m.  The  liner  feedback  configuration  used  in  the 
supercavitating  systems  is  again  utilized  here  (Eqs.  (2.6)).  Simulation  results  are 
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Figure  4.4:  Simulation  results  for  partial  cavity  model  with  a  =  0.066925. 

presented  in  Figure  4.4.  Since  there  is  no  planing  along  the  rear  of  the  vehicle,  it  is 
unable  to  support  itself  at  the  rear  (by  using  this  particular  control),  and  the  rear 
sinks  generating  a  high  pitch  angle  and  a  high  sideslip  angle  (as  shown  by  the  values 
for  the  transverse  velocity  w). 

Since  this  is  the  largest  cavity  that  can  be  converged  upon  by  using  the  partial 
cavity  numeric  model,  larger  partial  cavities  are  approximated  by  using  the  super- 
cavitating  numeric  model.  This  approach  may  not  produce  an  accurate  cavity  shape 
since  it  does  not  consider  the  vehicle  body  in  the  flow.  But,  cavities  generated  us¬ 
ing  this  approach  can  be  used  as  representative  partial  cavity  shapes  to  illustrate 
potential  partially  cavitating  vehicle  dynamics.  The  supercavitating  numeric  cavity 
model  is  run  at  low  speeds  (compared  to  the  supercavitating  vehicle  simulations), 
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Figure  4.5:  Simulation  results  obtained  by  using  numeric  supercavitatiug  model  to 
estimate  partial  cavity  shape  at  a  =  0.043. 

which  will  mean  the  generation  of  generate  small  cavities.  The  small  cavities  are 
then  truncated  to  only  sections  where  the  cavity  radius  is  larger  then  the  vehicle 
body.  The  truncated  supercavity,  is  then  treated  as  an  estimate  of  a  partial  cavity. 
Results  of  a  simulation  obtained  for  cr  =  0.043  is  shown  in  Figure  4.5.  This  corre¬ 
sponds  to  a  cavity  length  of  1.5106  m  (close  to  the  length  of  the  vehicle  at  L  =  1.8 
m).  Here,  limit-cycle  motion  is  observed. 

By  running  the  simulation  with  a  slightly  smaller  cavity,  at  a  =  0.046,  with 
1.3930  m  of  the  vehicle  unwetted,  the  rear  of  the  vehicle  is  again  unable  to  be 
supported  by  the  planing  (see  Figure  4.6). 
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Figure  4.6:  Simulation  results  obtained  by  using  numeric  supercavitatiug  model  to 
estimate  partial  cavity  shape  at  cr  =  0.046. 

4.1.4  Fin  Feedback 

For  the  partially  cavitating  vehicles,  the  planing  forces  are  not  rearward  enough 
to  support  the  vehicle.  A  simulation  run  with  a  passive  fin  input  of  Se  =  0.1  rad, 
is  considered  (with  the  same  the  linear  cavitator  feedback  control).  As  mentioned 
in  the  previous  chapter,  this  is  approximately  the  £n  angle  required  to  statically 
support  the  rear  of  a  supercavitating  vehicle.  The  simulation  results  are  presented 
in  Figure  4.7.  For  the  supercavitating  system,  the  vehicle  was  able  to  eventually 
stabilize  inside  the  cavity  using  this  passive  £n  input  (no  planing  support  required). 
However,  for  the  partially  cavitating  system,  even  with  the  buoyancy  forces  helping 
to  support  the  vehicle  (as  compared  to  the  supercavitating  case),  the  vehicle  is  still 
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Figure  4.7:  Simulation  results  obtained  for  the  partial  cavity  model  a  =  0.066925, 
with  passive  fin  input  of  5e  =  0.1  rad. 

unable  to  support  itself.  This  may  be  attributed  to  the  fact  that  for  the  partially 
cavitating  system,  the  planing  forces  and  moments  are  insufficient  to  resist  transient 
motions,  whereas  with  the  supercavitating  system,  planing  is  able  to  reject  high  pitch 
rates  and  vehicle  sideslip  angles  before  the  vehicle  stabilizes  inside  the  cavity. 

A  modified  linear  fin  feedback  in  addiition  to  the  linear  cavitator  feedback  is 
considered  next.  A  simulation  conducted  using  the  control  actuation  given  by  Eq. 
(4.10)  is  shown  in  Figure  4.8.  Similar  to  the  reasoning  used  for  the  choice  of  feedback 
states  in  Eqs.  (2.6),  the  fin  feedback  is  based  on  a  practically  measurable  state  (as 
opposed  to  vehicle  side  slip  w).  Here,  the  vehicle  clearly  stabilizes  with  very  small 
sideslip  angle  (moving  predominantly  forward  with  respect  to  the  vehicle  axis). 
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Figure  4.8:  Simulation  results  obtained  for  the  partial  cavity  model  with  a  = 
0.066925,  using  cavitator  and  fin  linear  feedback. 


Se  =  .12 +  .3  q 

Sc  =  l5z-306-.3q  (4.10) 

The  Matlab  code  used  for  the  partially  cavitating  vehicle  dynamics  model  is 
included  in  Appendix  B. 

4.1.5  Summary 

The  partial  cavitating  vehicle  dynamics  model  was  generated  to  explore  vehicle 
motions  present  during  this  type  of  operation.  A  numeric  partial  cavity  model  is 
used,  but  this  model  only  provides  solutions  for  limited  length  cavities.  Longer 
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partial  cavities  are  estimated  by  using  truncated  versions  of  supercavity  solutions. 
Vehicles  operating  under  partially  cavitating  conditions  are  modeled  here  by  adding 
buoyancy  forces  and  simplified  added  mass  expressions.  Because  the  cavity  closes 
further  forward  on  the  body,  compared  to  the  supercavitating  systems,  planing 
does  not  provide  sufficient  resistance  (for  smaller  partial  cavities)  to  high  vehicle 
sideslip  angles  without  the  aid  of  feedback  fin  control.  A  future  goal  would  be  the 
capability  of  simulating  full  vehicle  missions,  from  fully  wetted,  to  cavity  growth 
(partial  cavitation),  to  full  supercavitation  (and  potentially  transitioning  between 
the  different  operating  conditions).  Unfortunately  at  the  time  of  this  work,  there  are 
no  unsteady  partial  cavity  models.  The  work  presented  here  provides  an  initial  step 
towards  modeling  full  vehicle  missions  that  include  cavity  growth  (and  collapse). 

4.2  Delayed  Cavity  with  Non-Steady  Planing  Forces 

In  much  of  the  previous  research  (and  in  all  of  the  aforementioned  models),  the 
planing  force  modeling  is  based  on  the  assumption  that  steady  planing  is  neeeded 
to  simplify  the  planing  force  calculations.  With  this  assumption,  vehicle  motions 
into  or  out  of  the  cavity  are  ignored.  The  simplification  pertains  to  the  immersion 
rate  term  h.  For  steady  planing,  the  immersion  rate  at  any  section  of  planing  can 
be  represented  as  h  —  Vt  ■  sin(a),  where  a  is  the  angle  between  the  body  and  the 
cavity,  and  Vt  in  this  expression  represents  the  total  vehicle  speed  (Figure  4.9).  This 
immersion  rate  accounts  for  the  fact  that  the  body  is  moving  along  the  cavity  axis 
at  an  angle  but  does  not  include  any  radial  motions  of  the  body  into  the  fluid.  With 
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Figure  4.9:  Diagram  of  a  cylinder  planing  on  a  cylindrical  surface. 

the  steady  planing  assumption,  the  planing  force  becomes  only  a  function  of  the 
body’s  position  with  respect  to  the  cavity,  and  therefore  it  can  be  considered  as 
having  no  damping  relationship  in  terms  of  the  vehicle’s  motion. 

As  demonstrated  in  the  previous  models,  cavity  shape  and  location  predictions 
can  have  a  significant  effect  on  the  resulting  dynamics.  Two  methods  of  modeling 
cavity  position  and  orientation  (with  respect  to  the  vehicle  body)  are  shown  in 
Figure  4.10.  The  instantaneous  approach  is  utilized  in  all  of  the  previous  models. 
With  this  approach,  the  cavity  position  and  orientation  are  calculated  based  on  the 
current  cavitator  position  and  velocity  direction.  In  this  approach,  one  approximates 
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Instantaneous  approach:  Cavity  estimate  based 
on  current  cavitator  position  and  velocity  data 


cavitator  position  and  velocity  data 

Figure  4.10:  Two  methods  of  modeling  cavity  position  and  orientation. 

the  cavity  at  the  rear  of  the  vehicle  (where  planing  occurs)  as  a  cavity  that  would 
have  been  generated  if  the  cavitator  had  been  moving  along  it’s  current  velocity 
direction  up  to  it’s  current  position.  This  approximation  works  well  at  capturing 
the  dynamics  of  the  system  when  the  vehicle  speeds  are  high.  However,  the  cavity 
at  the  rear  of  the  vehicle  is  actually  generated  by  previous  motions  of  the  cavitator 
through  the  fluid  (which  may  not  coincide  with  the  current  position  and  orientation 
of  the  cavitator).  As  such,  a  more  appropriate  method  of  representing  the  cavity 
for  planing  is  to  model  the  cavity  centerline  and  orientation  based  on  a  previous 
cavitator  position  and  velocity  information.  This  is  shown  as  the  delayed  approach 
in  Figure  4.10. 

In  this  section,  a  model  is  presented  which  uses  both  a  delayed  cavity  approach 
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Cavity  Centerline 


Figure  4.11:  Cavity  and  vehicle  centerlines  for  the  delayed  case. 

to  generate  cavity  location  and  orientation  for  planing,  as  well  as  a  planing  force 
formulation  that  includes  vehicle  motions  into  and  out  of  the  fluid.  For  simplification 
and  computational  purposes,  only  cylindrical  cavities  are  considered.  However,  the 
approach  may  be  extended  to  account  for  non-cylindrical  cavity  planing. 

4.2.1  Description  of  Immersion  Terms 

A  depiction  of  the  cavity  and  vehicle  centerlines  for  an  arbitrary  vehicle  posi¬ 
tion  is  shown  in  Figure  4.11.  The  portion  of  the  cavity  that  interacts  with  the  vehicle 
is  generated  by  previous  positions  and  orientations  of  the  cavitator.  When  using  a 
cylindrical  cavity  approximation,  the  cavity  radius  and  location  can  be  calculated 
by  using  a  single  previous  cavitator  position  and  orientation  (a  single  delay).  The 
delay  is  taken  as  r  seconds,  which  corresponds  to  the  amount  of  time  it  took  for 
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the  cavitator  to  move  from  the  position  of  interest  (where  the  planing  is  occurring 
at  the  back  of  the  vehicle)  to  where  it  is  currently.  In  the  figure,  the  cavity  center 
is  located  at  the  previous  nose  position  (; xT,zT ).  The  cavity  expands  in  a  radial 
direction  perpendicular  to  the  velocity  direction  at  the  previous  time,  making  the 
cavity  axis  parallel  to  the  delayed  velocity  direction.  The  cavity  angle  with  respect 
to  horizontal  is  denoted  as  9cT  and  can  be  expressed  as  shown  in  Eq.  (4.11).  In  this 
expression,  0T  represents  the  delayed  body  orientation  with  respect  to  horizontal, 
and  tan~1(wT/V )  represents  the  delayed  velocity  orientation  with  respect  to  the 
body. 


0cT  =  9t  —  tan  1(wT/V)  (4.11) 

The  relative  position  of  the  rear  of  the  body  with  respect  to  the  cavity  cen¬ 
terline  is  a  function  of  both  a  translation  and  a  body  rotation.  The  translation  of 
the  current  position  of  the  nose  with  respect  to  the  cavity  centerline  is  expressed  as 
b.  The  relative  angle  between  the  current  body  centerline  and  the  cavity  centerline 
can  be  expressed  as  9(t)  —  9cT,  where  9(t)  is  the  current  body  orientation  with  re¬ 
spect  to  horizontal.  The  displacement  due  to  body  rotation,  c,  can  be  expressed  as 
c  =  a  ■  tan{9(t)  —  9cr ).  The  radial  displacement  of  the  body  centerline  at  the  rear 
with  respect  to  the  cavity  at  (xT,  zT ),  is  simply  b  +  c. 

The  immersion  depth  h  can  then  be  described  as  given  in  Eq.  (4.12),  where 
A  =  Rc  —  R ;  this  is  a  simplified  expression  for  immersion  only  along  one  surface. 
The  immersion  rate  can  then  be  generated  by  differentiating  Eq.  (4.12)  which  yields 
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Eq.  (4.13).  The  delay  terms  are  treated  as  having  no  dependence  on  time.  They 
relate  only  to  a  specific  instance  in  time  which  is  used  to  generate  an  instantaneous 
cavity  orientation  in  space  for  use  in  the  calculation  of  the  planing  forces.  With 
the  exception  of  expansion  and  contraction  along  the  cavity  radial  direction,  this 
instantaneous  cavity  does  not  move  or  change.  This  corresponds  to  the  physical 
understanding  that  the  cavity  is  not  moving  in  space  once  created;  that  is  it  is 
simply  expanding  or  contracting. 

h  —  a  -  tan(9(t)  —  dcT )  +  b  —  A  (4-12) 

h  —  a  ■  sec2(d(t )  —  dcT)  ■  q(t )  +  a  ■  tan(9(t )  —  9cr)  +  b  —  Rc  (4.13) 

The  terms  a  and  b  represent  the  motion  of  the  vehicle  nose  relative  to  the 
fixed  cavity.  These  terms  are  the  axial  and  radial  (with  respect  to  the  cavity  axis) 
components  of  the  nose  velocity.  Going  back  to  the  immersion  rate  expression  in 
Eq.  (4.13),  the  first  term  relates  to  the  rotation  of  the  body  into  the  cavity,  the 
second  term  relates  to  the  fact  that  the  body  is  moving  through  the  cavity  with  a 
relative  angle,  the  third  term  relates  to  the  rigid  body  motion  of  the  vehicle  into  the 
fluid,  and  the  last  term  relates  to  the  cavity  radial  growth  rate. 

The  parameters  a  and  b  can  be  solved  by  using  geometry.  In  Figure  4.12, 
the  orientation  of  a  and  b,  along  with  (xT,zT)  and  (x(t),  z(t)),  are  shown.  The  line 
segment  that  joins  (xT,  zT),  and  (x(t),  z(t)),  creates  an  angle  of  9X  =  tan -1 
with  the  horizontal;  the  term  is  zT—z(t)  since  z  is  positive  in  the  downward  direction. 
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Cavity  Centerline 


Figure  4.12:  Parameters  a  and  b  in  relation  to  cavity  centerline. 

The  relative  angle  that  this  segment  creates  with  the  cavity  axis  is  simply  9cr  —  6X. 
The  parameters  a  and  b  can  then  be  expressed  as  in  Eqs.  (4.14)-(4.15).  By  expanding 
the  sin  and  cos  terms,  the  expressions  can  be  simplified  to  Eqs.  (4.16)-(4.17).  The 
rate  of  change  can  then  be  expressed  as  given  in  Eqs.  (4.18)-(4.19).  As  described 
earlier,  these  terms  can  also  be  considered  as  the  axial  and  radial  components  (with 
respect  to  the  cavity  axis)  of  the  vehicle  velocity  at  the  nose. 


a  = 


b  = 


(zT  -  z(t ))2  +  ( x(t )  -  xT )2  •  cos  j 
\J (zT  -  z(t))2  +  (x(t)  -  xT)2  ■  sin  j 
a  =  cos{6cT)(x{t ) 
b  =  sin{6cT)(x[t) 


—  xT )  +  sin(9cT)(zT  —  z(t )) 


xT)  —  cos(9ct)(zt  —  z(t)) 


(4.14) 

(4.15) 

(4.16) 

(4.17) 


a  =  x  ■  cos[9ct )  —  z  ■  sin{9cT) 


(4.18) 
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b  —  x  ■  sin{6cT)  +  z  ■  cos(9cT ) 


(4.19) 


As  an  aside,  the  impacting  planing  force  expressions  can  be  solved  for  the  non- 
delayed  case.  In  the  case  with  no  delay,  the  cavity  is  directly  related  to  the  current 
conditions  at  the  nose,  and  the  axis  of  the  cavity  is  oriented  along  the  current 
velocity  direction.  For  no  delay,  0T  =  6{t),  wT  =  w(t),  zT  =  z(t),  and  b  =  b  =  0.  The 
expression  for  the  immersion  depth  becomes  hnon_deiay  =  a  ■  w/V  —  A  where  V  is 
the  forward  vehicle  speed  (as  used  in  the  dynamics  modeling).  The  immersion  rate 
simplifies  to  hnon _deiay  =  a  ■  sec2 {tan~l {w /V))  ■  q  +  a  ■  w/V  —  Rc,  where  a  is  total 
vehicle  speed. 

4.2.2  Integration  into  Dynamics  Model 

In  order  to  incorporate  the  delay,  the  overall  vehicle  path  needs  to  be  accurately 
tracked  in  the  inertial  frame.  The  small  angle  assumptions  can  be  removed  from 
the  propagation  of  the  depth  state  z  and  an  additional  state  for  the  x  position  can 
be  added  to  the  equations  of  motion  represented  by  Eqs.  (2.1)  and  (2.2).  The 
expressions  for  i  and  x  are  shown  as  follows. 


z  —  w  ■  cos (9 )  —  V  ■  sin(9 ) 

(4.20) 

x  —  V  ■  cos(9 )  +  w  ■  sin{6 ) 

(4.21) 

The  planing  force  is  calculated  from  the  Paryshev  representation  [37],  which 
is  solved  for  the  cylinder  on  cylinder  case.  The  resulting  planing  force  can  then  be 
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represented  as  shown  in  Eq.  (4.22),  where  A  =  Rc  —  r,  and  ho  is  the  immersion 
depth  at  the  aft  of  the  vehicle,  as  shown  in  Figure  4.9. 

1  /  A2  \ 

FP  =  7T  pr2h2 - 1  —  — - -  (4.22) 

tan(a)  y  (h0  +  A)2y 

4.2.3  Simulation  Results 

For  the  simulation  runs,  the  vehicle  parameters  were  chosen  to  match  those 
used  in  the  previous  runs  with  m  =  2,  Rn  =  0.0191  m,  R  =  0.0508  m,  L  =  1.8  m, 
n  =  0.5,  and  Cx o  =  0.82.  Feedback  control  is  again  described  by  Eqs.  (2.6)  where 
the  fins  are  assumed  to  be  passive,  while  the  cavitator  utilizes  linear  state  feedback. 
The  nominal  value  for  the  delay  simulations  is  chosen  as  t  =  L/V,  the  approximate 
amount  of  time  it  takes  for  the  vehicle  nose  to  travel  one  body  length.  For  low 
cavitation  numbers  (high  speeds),  the  system  shows  a  stable  equilibrium  response. 
Similar  to  the  other  dynamics  models  with  the  same  feedback  formulation,  this 
system  exhibits  oscillations  as  the  cavitation  number  is  increased  (speed  is  decreased) 
and  the  cavity  to  body  clearance  tightens. 

A  simulation  run  conducted  for  a  =  0.0241  (V  =  83.68  m/s )  is  shown  for  the 
delayed  system  in  Figure  4.13.  This  particular  value  of  a  is  chosen  to  be  just  below 
the  critical  cavitation  number  where  limit-cycle  motion  is  observed. 

A  simulation  run  for  the  non-delayed  system  at  the  same  cavitation  number 
a  =  0.0241  is  shown  in  Figure  4.14.  ffere  the  non-delayed  system  exhibits  limit 
cycles.  It  should  be  noted  that  for  high  speeds  (low  cavitation  numbers),  the  non- 
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Figure  4.13:  System  response  results  for  cr  =  0.0241  and  delayed  system. 

delayed  system  does  demonstrate  stable  behavior.  So  when  starting  from  sufficiently 
high  speeds,  as  the  cavitation  number  is  increased,  both  the  delayed  and  non-delayed 
system  transition  from  stable  to  limit-cycle  motions.  However,  the  transition  for 
the  non-delayed  system  occurs  slightly  earlier.  It  is  within  this  window  of  cavitation 
numbers,  where  the  equilibrium  position  of  the  delayed  system  is  stable  and  the 
non-delayed  system  is  unstable,  that  the  delay  can  be  considered  stabilizing.  This 
is  in  contrast  to  previous  findings  where  the  delay  can  be  shown  to  destabilize  the 
system 

If  the  delay  is  taken  as  a  parameter,  and  the  forward  velocity  is  held  constant 
at  V  —  83.68  m/s,  the  steady-state  behavior  of  the  system  follows  the  plots  shown 
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Figure  4.14:  System  response  results  for  a  =  0.0241  and  non-delayed  system. 


in  Figure  4.15.  Here,  the  delay  is  varied  from  close  to  0  s  (no  delay)  to  0.01936 
s  (approximately  90%  of  the  nominal  delay  value  of  L/V).  The  system  can  be 
shown  to  transition  from  limit-cycle  motion  to  asymptotic  stability  as  the  delay  is 
increased.  By  varying  the  delay  with  a  fine  timestep  of  less  then  0.00001  s,  the 
critical  value  of  r  =  0.01930  s  can  be  found,  where  limit  cycles  exist  for  r  <  0.01930 
s. 
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Figure  4.15:  Effect  of  delay  on  steady-state  system  response. 
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Chapter  5 


Vehicle  Maneuvering  Using  Numeric  Optimal  Control  Approach 

In  this  chapter,  a  framework  is  provided  for  analyzing  maneuvering  of  non¬ 
smooth  vehicle  systems.  Maneuverability,  for  nonlinear  and  particularly  for  non¬ 
smooth  vehicle  systems,  can  be  difficult  to  characterize  since  “all  out”  or  fully 
saturated  control  inputs  do  not  necessarily  define  the  envelope  of  vehicle  motion 
capabilities.  For  example,  in  the  case  of  supercavitating  vehicles,  a  fully  saturated 
control  input  will  push  the  system  into  tailslap  which  negatively  affects  the  position 
control  of  the  vehicle.  In  this  work,  a  numeric,  direct  optimal  control  approach  has 
been  used  for  generating  optimal  control  inputs  for  defined  vehicle  maneuvers.  This 
approach  is  applied  to  the  supercavitating  vehicle  systems  presented  in  the  previ¬ 
ous  chapters.  The  supercavitating  systems  are  complex  examples  of  non-smooth 
dynamic  systems,  they  serve  as  good  candidates  for  numerical  approaches. 

5.1  General  Approach 

Due  to  the  complexity,  and  non- differentiable  nature  of  non-smooth  systems, 
the  direct  method  for  optimal  control  was  chosen.  In  this  method,  controller  inputs 
are  determined  numerically.  The  control  inputs  for  discrete  segments  of  time,  or 
parameters  that  are  used  to  define  a  control  input  function,  are  treated  as  variables 
in  an  optimization  scheme.  The  system  dynamics  can  be  directly  integrated  by 
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Figure  5.1:  Diagram  of  discretization  method  for  the  optimization  strategy  used  to 
generate  controller  input. 

using  the  governing  equations  of  motion.  The  maneuver  itself  can  be  treated  as  a 
constraint  on  the  optimization  formulation,  constraining  the  dynamics  to  perform  a 
desired  task.  The  control  input  (s)  can  then  be  optimized  for  total  time  (which  can 
also  be  treated  as  a  variable),  or  ending  state,  depending  on  the  type  of  maneuver 
considered. 

A  simple  example  of  an  application  would  be  a  maneuver  subject  to  optimiza¬ 
tion  for  total  time  of  maneuver  T,  which  has  end  point  constraints  (such  as  final 
position  and  orientation  constraints).  The  maneuver  can  then  be  discretized  into 
s  equal  length  time  segments  where  the  control  input  is  constant  over  each  time 
segment,  as  shown  in  Figure  5.1.  The  values  for  the  control  tq  for  %  —  1 . . .  s  as 
well  as  the  final  time  T  are  considered  as  variables  in  the  optimization  scheme.  The 
dynamics  can  then  be  directly  integrated  from  a  given  initial  state  Xq  and  with  the 
final  position  xs,  subject  to  the  constraints  determined  by  the  type  of  maneuver 
being  considered.  This  type  of  maneuver  is  illustrated  in  Figure  5.2. 

The  optimization  formulation  can  be  described  as  follows.  The  objective  func- 
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Figure  5.2:  Diagram  of  trajectory  generated  by  using  the  constant  control  inputs 
over  discretized  time  segments. 

tion  can  be  expressed  as  given  in  Eq.  (5.1).  This  function  is  subject  to  the  constraint 
shown  in  Eq.  (5.2),  where  s  is  the  number  of  time  intervals  and  Xf  is  defined  by  the 
maneuver.  The  state  values  at  end  of  each  of  the  time  intervals  can  be  expressed  as 
shown  in  Eq.  (5.3),  for  i  =  l...s  and  Xq  a  given.  Additional  constraints  can  also  be 
applied  to  bound  the  control  inputs  and/or  states. 


THXTlu  1  ? u2,---,us  ,t(T) 

(5.1) 

Xs  =  Xf 

(5.2) 

rt=iT/s 

Xi  —  1  F[t,x,Ui)dt 

Jt=(i—l)T/s,Xi-i 

(5.3) 

An  important  consideration  in  this  formulation  is  that  this  formulation  in- 
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Figure  5.3:  Simple  point  mass  controlled  system. 

herently  only  considers  trajectories  that  the  system  is  capable  of  achieving  (since 
the  dynamics  are  directly  integrated).  This  eliminates  the  need  to  determine  a  set 
of  achievable  trajectories  as  done  with  some  path  planing  methods.  However,  end 
conditions  must  be  chosen  with  care,  since  an  un-achievable  end  condition  means 
that  there  is  an  empty  feasible  set  of  variables  for  the  optimization  formulation. 


5.1.1  Application  to  Simple  System 

To  better  illustrate  the  optimal  control  method,  it  is  applied  to  a  planar  point 
mass  system.  The  system  is  shown  in  Figure  5.3,  wherein  the  inputs  are  Fx  and  Fy. 
They  represent  the  forces  along  the  x  and  y  directions  respectively.  The  equations 
of  motion  are  given  by  Eq.  (5.4),  where  the  state  vector  is  given  by  (x,y,x,y). 
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t  (s) 

x  (m) 

y  (■ m ) 

x  (m/s) 

y  (m/s) 

Fx  (N) 

Fy  (A0 

0.5 

0.5 

0.5 

2.2 

2.2 

5.0 

5.0 

0.9 

2.0 

2.0 

4.5 

4.5 

5.0 

5.0 

Table  5.1:  Optimal  solution  for  point  mass  system  for  starting  point 

(0.0m,  0.0m,  0.0 m/s,  0.0 m/s)  and  end  condition  (2.0m,  2.0m,  _,  _)  and  s  —  2. 

The  control  forces  are  bounded  to  [— 5N,  hN].  The  control  inputs  are  calcu¬ 
lated  by  using  the  constrained  optimizer  in  Matlab,  fmincon.  Simple  maneuvers  are 
considered  first,  with  s  —  2  and  an  initial  condition  of  xq  =  (0.0m,  0.0m,  0.0 m/s,  0.0 m/s) 
(which  corresponds  to  an  initial  position  at  the  origin,  with  zero  speed).  A  so¬ 
lution  for  a  maneuver  that  is  only  defined  by  a  final  position  is  shown  in  Ta¬ 
ble  5.1,  where  ( Xf,yf )  =  (2.0m,  2. 0m)  (any  velocity  at  final  position  allowed). 
The  optimal  solution  for  a  fastest  time  maneuver  is  a  “full-on”,  u  =  5.0,  style 
control.  A  second  maneuver  that  defines  both  the  final  position  and  final  ve¬ 
locity  is  shown  in  Table  5.2.  Here,  the  maneuver  is  defined  with  the  end  point 
(xf,  i/f,  Xf,  y'f)  =  (2.0  m,  2.0  m,  O.Om/s,  0.0m/s),  which  can  be  thought  of  as  a  move- 
to  and  stop.  The  optimal  control  inputs  for  this  maneuver  are  a  “full-on”  ( u  =  5.0) 
time  segment,  followed  by  a  “full-stop”  ( u  =  —5.0)  control.  The  solution  for  these 
symmetric  position  and  velocity  end  conditions  are  not  surprising,  and  they  can 
easily  be  determined  by  inspection.  The  optimal  control  process  discussed  here  can 
be  applied  to  more  complicated  maneuvers. 

The  optimal  time  solution  for  a  maneuver  to  the  end  point  (xf,yf,Xf,yf)  = 
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t  (s) 

x  ( m ) 

y  (■ m ) 

x  (m/s) 

y  (m/s) 

Fx  (N) 

Fy  (A0 

0.3 

0.3 

0.3 

1.6 

1.6 

5.0 

5.0 

0.6 

1.0 

1.0 

3.1 

3.1 

5.0 

5.0 

.09 

1.8 

1.8 

1.6 

1.6 

-5.0 

-5.0 

1.3 

2.0 

2.0 

0.0 

0.0 

-5.0 

-5.0 

Table  5.2:  Optimal  solution  for  point  mass  system  for  startint  point 

(0.0m,  0.0m,  O.Om/s,  O.Om/s)  and  end  condition  (2.0m,  2.0m,  O.Om/s,  O.Om/s)  and 
s  =  2. 

(1.0m,  7.0m,  — l.Om/s,  2.0m/s)  is  less  intuitive,  and  the  best  result  found  from  the 
optimization  process  is  shown  in  Table  5.3.  The  trajectory  as  well  as  the  time 
histories  of  the  sates  are  also  shown  in  Figure  5.4.  The  resulting  motion  can  be  seen 
to  move  around  in  the  x  direction;  this  is  due  to  the  fact  that  the  limiting  factor 
for  fastest  time  maneuver  was  the  motion  in  the  y  direction.  Although  the  solution 
here  is  more  difficult  to  ascertain  via  inspection,  simple  systems  are  only  useful  to 
validate  the  numeric  optimal  control  approach.  This  type  of  approach  is  intended 
for  more  complex  systems. 

5.2  Maneuvering  with  Cylindrical  Dive-Plane  Models 

The  optimal  control  approach  is  applied  to  a  supercavitating  vehicle  system. 
The  first  system  considered  is  based  on  the  original  cylindrical  cavity,  dive-plane 
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Time  (s) 


Time  (s) 


Figure  5.4:  Best  solution  for  point  mass  system  for  starting  point 

(0.0m,  0.0m,  0.0 m/s,  0.0 m/s)  and  end  condition  (1.0m,  7.0m,  —1.0 m/s,  2.0m/s)  and 
s  =  10. 
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step  : 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

F,  (TV) 

5.0 

4.5 

-3.4 

-5.0 

-5.0 

5.0 

1.1 

2.8 

-5.0 

-5.0 

F,  (AO 

5.0 

5.0 

5.0 

5.0 

5.0 

4.8 

-5.0 

-5.0 

-5.0 

-5.0 

T  =  2.0 

Table  5.3:  Best  solution  for  point  mass  system  for  starting  point 

(0.0m,  0.0m,  O.Om/s,  O.Om/s)  and  end  condition  (1.0m,  7.0m,  — l.Om/s,  2.0m/s) 
and  s  =  10. 

model,  which  utilizes  the  Hassan  planing  force  model  (as  used  in  reference  [8]). 
This  model  was  chosen  to  implement  first  because  of  it’s  simplicity  (compared  to 
the  subsequent  modified  models),  while  the  system  still  exhibited  behavior  such  as 
limit-cycle  motions,  and  unstable  responses. 

Similar  to  the  delayed  system,  maneuvering  considerations  require  accurate 
tracking  of  the  inertial  position.  An  additional  state  for  the  x  coordinate  is  added, 
and  the  expression  for  the  propagation  of  z  is  adjusted  to  remove  the  small  angle 
assumptions.  The  state  space  representation  of  the  inertial  states,  x  and  z,  are 
shown  in  Eqs.  (5.5)-(5.6). 


x  —  V  ■  cos{6 )  +  w  ■  sin{6 ) 

(5.5) 

z  —  w  ■  cos (9 )  —  V  ■  sin{6 ) 

(5.6) 

Additional  constraints  are  applied  to  the  supercavitating  vehicle  system.  Ac¬ 
tuation  angles  for  both  the  fin  and  cavitator  control  surfaces  are  bounded,  and 
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these  can  be  applied  as  direct  bounds  on  the  optimization  variables.  Additionally, 
since  the  planing  force  model  is  only  valid  for  small  immersions,  state  limits  on  the 
transverse  speed,  |w|  <  6  m/s,  are  applied.  The  constraint  equations  can  then  be 
specified  as  shown  in  Eq.  (5.7)  and  Eq.  (5.8);  ceq  =  0  constrains  the  end  condition, 
and  q  <  0  is  used  to  apply  bounds  on  the  transverse  velocity  at  the  end  of  each  time 
interval.  Unless  otherwise  specified,  for  these  runs,  a  =  0.03,  which  corresponds  to 
a  velocity  (cavity  size)  where  the  system  experiences  limit-cycle  motions  and  un¬ 
stable  behavior.  All  runs  start  with  an  initial  condition  of  straight  and  level  flight 
originating  from  (z0,  xq,  wq,  qo)  =  (0.0/n,  0.0m,  O.Om/s,  O.Orad/s).  The  optimization 
(unless  otherwise  specified)  was  carried  out  by  using  an  off  the  shelf  constrained 
optimizer,  fmincon,  in  Matlab  (and  specifically  the  interior  point  algorithm  for  this 
function). 


‘eq  *£/  %s 

(5.7) 

Ci  =  \wi\  -  6 

Vi  =  1 ... s 

(5.8) 

If  the  optimal  control  approach  is  applied  without  the  aid  of  feedback  control, 
certain  solutions  can  be  found.  An  example  of  a  slight  dive  maneuver  is  shown  in 
Figure  5.5.  This  is  a  maneuver  to  Zf  =  2.0  m  and  Wf  =  0.0  m/s,  which  corresponds 
to  a  dive-to-depth,  and  an  alignment  of  the  body  with  the  total  velocity  direction. 
For  this  run  s  =  8  and  an  initial  guess  of  zero  cavitator  input  and  zero  fin  input  is 
applied.  In  this  situation,  the  optimizer  was  capable  of  finding  a  solution  with  a  fea¬ 
sible  end  condition.  However,  if  the  discretization  was  increased  to  s  —  16,  a  feasible 


solution  by  using  the  same  starting  condition  could  not  be  found.  This  finding  was 
surprising,  since  the  solution  for  s  =  8  is  clearly  feasible  for  s  =  16.  Additionally, 
feasible  solutions  for  other  similarly  simple  maneuvers  (including  straight  and  level 
flight)  conld  not  be  found,  even  with  good  initial  conditions.  The  issue  relates  to 
the  fact  that  the  system  without  feedback  control  is  inherently  unstable  (as  shown 
in  the  previous  work  and  in  Chapter  2).  As  the  optimizer  searches  for  solutions  that 
satisfy  the  end  conditions,  the  system  can  easily  transition  into  instability.  When 
no  feedback  control  is  utilized,  the  optimal  control  method  is  tasked  with  generat¬ 
ing  a  controller  that  stabilizes  the  system,  while  maneuvering  it  to  the  proper  end 
condition.  However,  with  the  exception  of  some  specific  cases,  the  coarseness  of  the 
discretization  used  for  the  optimal  control  approach  was  found  to  be  incapable  of 
rejecting  planing  instabilities. 

5.2.1  Inner-Loop  and  Outer-Loop  Control  Schemes 

A  feedback  controller  can  be  added  to  the  system.  The  feedback  controller  can 
be  considered  as  an  “inner-loop”  controller  that  is  used  to  help  reject  fast  timescale 
instabilities  (such  as  with  the  planing  force),  while  the  control  inputs  generated  by 
the  optimization  can  be  considered  as  “outer-loop”  control  that  guides  the  vehicle 
through  the  desired  maneuver  to  the  proper  end  condition.  A  diagram  of  how 
the  two  controllers  are  integrated  is  shown  in  Figure  5.6.  Here,  the  optimization 
formulation  plays  the  role  of  the  motion  planner. 

By  utilizing  the  liner  feedback  law  shown  in  Eqs.  (2.6),  the  system  is  stabilized 


87 


Figure  5.5:  Trajectory  for  dive  maneuver  to  Zf  =  2.0  m  and  Wf  =  0.0  m/s  with  no 
feedback  control  and  s  =  8. 


to  stable  limit  cycles.  The  response  time  histories  for  a  solution  using  this  feedback 
law  for  straight  and  level  flight  is  shown  in  Figure  5.7,  along  with  the  determined 
controller  inputs  in  Figure  5.8.  ffere,  the  maneuver  requires  an  end  condition  of 
(xf,Zf,Wf)  =  (40.0m,  0.0m,  0.0m/ s)  with  all  other  final  states  unconstrained.  By 
using  an  initial  guess  of  constant  fin  input  of  Sf  =  0.1  rad  (approximately  what  is 


Outer  loop 


Figure  5.6:  Depiction  of  inner-loop  and  outer-loop  controllers. 


Figure  5.7:  Time  histories  for  maneuver  to  ( Xf,Zf,Wf )  =  (40.0m,  0.0m,  0.0m/ s) 
with  feedback  control  according  to  Eqs.  (2.6)  and  s  =  8. 

required  to  support  the  rear  of  the  vehicle),  the  optimization  scheme  was  easily  able 
to  find  a  good  set  of  control  inputs. 

One  issue  with  this  particular  form  of  feedback  is  that  it  directly  involves  the 
inertial  states  z  and  6.  Unfortunately,  this  greatly  limits  the  ability  to  perform 
maneuvers  aside  from  straight  and  level  flight.  Instead,  the  inertial  terms  can  be 
dropped  and  an  inner-loop  control  of  the  form  as  shown  in  Eq.  (5.9)  can  be  utilized. 
This  type  of  controller  still  aids  in  stabilizing  the  system  (using  measurable  states), 
while  allowing  for  motions  in  space. 
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Figure  5.8:  Outer-loop  control  inputs  for  maneuver  to  (. Xf,Zf,Wf )  = 

(40.0m,  0.0m,  0.0m/s)  with  feedback  control  according  to  Eqs.  (2.6)  and  s  =  8. 


By  using  kinner  =  —0.9,  a  dive  maneuver  to  Zf  =  2  m  and  Wf  =  0  m/s 
is  considered  (all  other  final  states  being  unconstrained).  The  resulting  state  and 
outer-loop  control  histories  for  the  best  solution  are  shown  in  Figures  5.9  and  5.10. 
An  initial  guess  of  zero  cavitator  or  hn  control  angles  is  utilized  here  with  s  =  8. 
Also,  unlike  the  case  with  no  inner-loop  control,  control  inputs  for  similar  simple 
maneuvers  can  be  easily  solved.  Deeper  dive  maneuvers  can  be  solved  relatively 
quickly  by  seeding  the  optimization  scheme  with  previous  solutions  for  shallower 
maneuvers.  Trajectories  for  dive  maneuvers  to  Zf  =  4.0,  8.0  and  20.0  m,  are  shown 
in  Figure  5.11  (again  with  Wf  =  0). 

With  the  inner-loop  stabilization,  more  complicated  maneuvers  can  also  be 
considered.  A  dive  followed  by  a  level-off  maneuver  can  be  described  with  Zf  =  20.0 
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Figure  5.9:  Time  histories  for  maneuver  to  ( Zf,Wf )  =  (2.0m,  0.0m/ s)  with  modified 
feedback  control  according  to  Eq.  (5.9),  kinner  =  —0.9,  and  s  =  8. 

m  and  Of  =  0.0  rad.  This  maneuver  is  solved  with  s  =  14  and  kinnner  =  —0.9,  and 
the  resulting  trajectory  and  outer- loop  control  solution  are  shown  in  Figures  5.12 
and  5.13,  respectively. 

5.2.2  Homing  maneuvers 

Move-to-point  maneuvers  are  also  of  interest  for  maneuvering.  These  types  of 
maneuvers  can  be  characterized  by  specifying  the  final  position  ( Zf,Xf ).  By  using 
the  results  from  the  previous  dive  maneuvers,  an  obtainable  move  to  point  marker 
was  chosen  at  ( Zf,Xf )  =  (20.0m,  80.0m).  The  resulting  control  solutions  with  s  =  14 
and  kinner  =  —0.7  are  shown  in  Figures  5.14  and  5.15.  Included  in  the  trajectory  plot 
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Figure  5.10:  Outer-loop  control  inputs  for  maneuver  to  ( Zf,Wf )  =  (2.0m,  0.0m/ s) 
with  modified  feedback  control  according  to  Eq.  (5.9),  kinner  =  —0.9,  and  s  =  8. 


Figure  5.11:  Trajectories  for  dive  maneuvers  to  (zf,Wf)  =  (4.0/n,  0.0m/s), 
(zf,Wf)  =  (8.0m,  0.0?n/s),  and  (zf,Wf)  =  (20.0m,  0.0/n/s),  with  modified  feedback 
control  according  to  Eq.  (5.9),  kinner  =  —0.9,  and  s  =  8. 
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Figure  5.12:  Trajectory  for  (zf,  Of)  =  (20.0m,  O.Orad)  with  modified  feedback  control 
according  to  Eq.  (5.9),  kinner  =  —0.9,  and  s  =  14. 


Time  (s) 


Figure  5.13:  Outer-loop  control  inputs  for  maneuver  to  (zf,0f)  =  (20.0m,  0. Or  ad) 
with  modified  feedback  control  according  to  Eq.  (5.9),  kinner  =  —0.9,  and  s  =  14. 
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of  Figure  5.14,  is  a  trajectory  showing  the  vehicle  path  with  inner-loop  control  only, 
as  well  as  vehicle  orientation  and  cylindrical  cavity  plots  at  specific  points  along  the 
actual  trajectory  (exaggerated  in  size  to  show  detail).  It  is  clear  that  the  vehicle  is 
planing  during  the  maneuver.  The  control  history  plot  in  Figure  5.15  also  includes 
a  total  cavitator  actuation  angle  plot  showing  the  combination  of  the  inner-loop 
and  outer-loop  control.  Similar  move  to  point  maneuvers  can  also  be  solved.  The 
resulting  best  solved  trajectory  for  a  climb  maneuver  to  (xf,  Zf)  =  (80.0m,  —20.0m) 
is  shown  in  Figure  5.16 


Figure  5.14:  Trajectory  for  ( Zf,Xf )  =  (20.0m,  80.0m)  with  modified  feedback  control 
according  to  Eq.  (5.9),  kinner  =  —0.7,  and  s  =  14. 
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Figure  5.15:  Control  inputs  for  maneuver  to  ( Zf,Xf )  =  (20.0m,  80.0m)  with  modified 
feedback  control  according  to  Eq.  (5.9),  kinner  =  —0.7,  and  s  =  14. 

Homing  maneuvers  with  a  circular  obstacle  or  “no-fly”  area  are  considered 
next.  If  the  obstacle  is  centered  at  ( x0bst ,  zobst )  with  radius  robst,  the  obstacle  can  be 
modeled  as  a  single  constraint  that  specifies  a  minimum  distance  from  the  obstacle 
center  to  the  closest  point  along  the  trajectory.  The  constraint  relationship  can 
be  expressed  as  in  Eq.  (5.10);  the  right-hand  side  expresses  the  closest  distance 
to  the  obstacle  center  over  the  entire  trajectory  (determined  by  the  integration  of 
the  equations  of  motion),  and  the  radius  of  the  obstacle  robst  specifies  the  minimum 
allowable  distance. 
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Figure  5.16:  Trajectory  for  ( Zf,Xf )  =  (— 20. 0m,  80.0m)  with  modified  feedback 
control  according  to  Eq.  (5.9),  kinner  =  —0.7,  and  s  =  14. 


robst  >  min  sj (x(t)  -  xobst )2  +  (z(t)  -  zobst )2  (5.10) 

A  climb  maneuver  to  (zf,Xf)  =  (—20.0m,  80.0m)  with  an  obstacle  is  centered 
at  (zobst,  x0bst)  =  (— 5.0m,  40.0m)  with  a  radius,  r0bst  =  10.0  m  is  then  considered. 
The  position  is  chosen  so  that  the  obstacle  is  located  in  the  best  found  move-to- 
point  trajectory,  which  has  been  generated  for  the  obstacle-free  case  (as  shown  in 
Figure  5.16).  The  resulting  trajectory  is  shown  in  Figure  5.17.  The  inner-loop  only 
trajectory  is  also  shown  along  with  the  vehicle-cavity  orientation  plots  (exaggerated 
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in  size).  As  expected,  the  best  discovered  trajectory  has  the  vehicle  operating  un¬ 
der  extreme  maneuvering  conditions  (as  shown  by  the  high  cavity  immersion),  and 
positions  the  path  very  close  to  the  obstacle  boundary.  The  response  histories  are 
plotted  in  Figure  5.18,  and  the  control  histories  are  shown  in  Figure  5.19. 


Figure  5.17:  Trajectory  for  ( Zf,Xf )  =  (—20.0m,  80. 0m),  obstacle  at  (zobst,xobst)  = 
(— 5.0m, 40.0m),  r0bst  =  10.0  m.  With  modified  feedback  control  according  to  Eq. 
(5.9),  kinner  0.7,  and  s  14. 

Additional  obstacles  can  be  also  be  accommodated  by  using  additional  con¬ 
straints.  A  run  to  a  further  point  of  ( Zf,Xf )  =  (— 50.0m,  200.0m)  is  shown,  with  two 
obstacles  at  (zobst,xobst)  1  =  (-5.0m,  40.0m),  and  (zobst,  xobst)2  =  (-10. 0m,  175.0m), 
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Figure  5.18:  Time  histories  for  ( Zf,Xf )  =  (—20.0m,  80.0m),  obstacle  at 
(zobsti  xobst)  =  (—5.0m,  40.0m),  r0bst  =  10.0  m.  With  modified  feedback  control 
according  to  Eq.  (5.9),  kinner  =  —0.7,  and  s  =  14. 


both  with  radius  rabst  =  15.0  m,  is  shown  in  Figure  5.20.  The  first  obstacle  was 
chosen  to  be  in  the  path  of  the  obstacle-free  best  discovered  solution,  and  the  second 
obstacle  was  chosen  to  be  in  the  path  of  the  single  obstacle  best  discovered  solution. 

Moving  end  points  can  also  be  considered,  and  they  are  of  interest  since  the 
vehicle  motion  does  not  occur  instantaneously.  A  moving  end  point  can  be  added 
by  using  an  expression  for  the  final  end  point  constraint  in  terms  of  the  final  time 
T.  An  example  is  generated  by  setting  (x/,  Zf)  =  (80  +  10T,  —20  —  10T)  m,  letting 
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Figure  5.19:  Outer-loop  control  inputs  for  ( Zf,Xf )  =  (—20.0m,  80.0m),  obstacle  at 
( ZobstiXobst )  =  (—5.0m, 40.0m),  r0bst  =  10.0  m.  With  modified  feedback  control 
according  to  Eq.  (5.9),  kinner  =  —0.7,  and  s  =  14. 


the  desired  end  point  start  at  ( Xf,Zf )  =  (80.0m,  —20.0m),  and  giving  a  velocity  of 
10  m/s  in  both  the  x  and  z  directions,  away  from  the  initial  position  of  the  vehicle. 
The  plot  of  the  resulting  trajectory  is  shown  in  Figure  5.21. 


5.3  Maneuvering  with  Non-Cylindrical  Dive-Plane  Models 

The  optimal  control  method  was  then  applied  to  the  non-cylindrical  dive-plane 
models.  As  presented  in  Chapter  3,  the  non-smooth  dynamics  is  complicated  due  to 
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Figure  5.20:  Trajectory  for  ( Zf,Xf )  =  (—50.0m,  200.0m),  obstacles  at  (z0bst,x0bst)  i  = 
(— 5.0m, 40.0m),  and  (z0bst,x0bst) 2  =  (—10.0m,  175.0m),  both  with  r0bst  =  15.0  m. 
With  modified  feedback  control  according  to  Eq.  (5.9),  kinner  =  —0.7,  and  s  =  34. 

non-constant  cavity  boundaries  (due  to  the  cavity  shift)  and  a  varying  planing  force 
function  (due  to  the  cavity  shape  change).  A  homing  maneuver  is  considered  first 
by  using  the  splined  cavity  model  (to  ensure  accurate  planing  area  predictions).  A 
solution  for  dive  maneuver  to  Xf  —  20.0  m,  Zf  =  80.0  m  with  0  =  0.0335,  is  shown 
in  Figure  5.22.  The  vehicle  orientation  and  non-cylindrical  cavity  plots  at  specific 
points  along  the  trajectory  are  also  shown,  and  are  exaggerated  in  size  to  show 
detail.  With  the  non-splincd  cavity  model,  the  computational  time  was  significantly 
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Figure  5.21:  Trajectory  for  ( Xf,Zf )  =  (80  +  10T,  — 20  —  10T)  m  with  modified 
feedback  control  according  to  Eq.  (5.9),  kinner  =  —0.7,  and  s  =  14. 

reduced,  but  larger  tolerances  on  the  end  condition  were  required. 

Incorporating  aggressive  maneuvering  with  the  non-cylindrical  models  be¬ 
comes  more  difficult.  Several  attempts  at  solving  the  single  obstacle  case  (similar  to 
the  cylindrical  case  shown  in  Figure  5.17)  yielded  no  feasible  solutions.  One  example 
where  the  constrained  optimizer  failed  to  find  a  feasible  solution  is  shown  in  Figure 
5.23.  Although  both  cylindrical  and  non-cylindrical  models  produced  similar  bifur¬ 
cation  behavior;  however,  the  parameter  values  at  which  the  bifurcations  occur  are 
different,  which  may  lead  to  significantly  different  maneuvering  capabilities  when 
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Figure  5.22:  Non-cylindrical  planing  with  splined  cavity  model,  trajectory  for 
(. Xf,Zf )  =  (80.0m,  —20.0m)  with  feedback  control  according  to  Eq.  (5.9),  a  = 
0.0335,  kinner  =  —0.9,  and  s  =  14. 

run  at  similar  parameter  values. 

5.3.1  Optimization  using  Penalty  Methods  and  Simple  Search  Algo¬ 
rithms 

In  order  to  investigate  the  maneuvering  capability  difference  between  the  non- 
cylindrical  and  cylindrical  models,  maximum  turn  maneuvers  are  considered.  A 
maneuver  that  maximizes  (or  minimizes)  Of  over  a  set  time  T,  originating  from 
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Figure  5.23:  Non-cylindrical  planing  with  splined  cavity  model;  failed  run  for 
(. Zf,Xf )  =  (—20.0m,  80.0m),  obstacle  at  (z0j,st,  xobst)  =  (—5.0m,  40.0m),  and  robst  = 

10.0  m. 

some  initial  condition  (such  as  (20,  %o,  u>o,  Qo)  —  (0.0m,  0.0m,  O.Om/s,  0.0rad/s))  can 
be  considered.  The  built  in  Matlab  function  fmincon  is  capable  of  solving  this  type 
of  maneuver  for  the  cylindrical  case;  however,  this  direct  constrained  approach  has 
difficulty  improving  upon  initial  guesses  for  the  non-cylindrical  models.  For  the 
non-cylindrical  planing  system,  the  maneuver  was  solved  by  using  a  non-gradient 
based  optimization  algorithm  with  the  function  fminsearch.  Although  these  types  of 
maneuvers  are  only  considering  one  specific  capability  measure,  the  differences  be¬ 
tween  the  results  for  the  non-cylindrical  and  cylindrical  models  were  not  significant. 
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This  finding  may  mean  that  the  incapability  of  the  constrained  optimizer  to  find  a 
feasible  solution  for  the  single  obstacle  maneuver  with  non-cylindrical  planing,  may 
be  more  related  to  the  complexity  of  the  feasible  domain  rather  then  the  incapability 
of  the  vehicle  to  maneuver. 

The  family  of  control  inputs  that  generate  trajectories  that  satisfy  the  maneu¬ 
ver  conditions  define  the  feasible  domain,  and  this  domain  may  include  unconnected 
sets.  Optimizing  in  this  domain  can  be  difficult  and  the  complexity  of  the  feasible 
domain  increases  as  the  dynamics  get  complex.  The  direct  constrained  optimization 
can  instead  be  replaced  by  a  penalty  type  objective  function.  Here,  the  constraints, 
rather  then  being  directly  enforced,  are  treated  as  penalties  on  the  objective  function 
in  an  unconstrained  optimization  scheme.  The  optimizer  progresses  simultaneously 
towards  feasible  and  optimal  solutions. 

For  the  example  of  the  supercavitating  vehicle  system,  the  objective  function 
can  be  re-written  as  Eq.  (5.11).  The  values  ct  represent  the  inequality  constraints 
(as  used  for  the  bounds  on  the  states  and  control),  and  the  values  ceqj  represent 
the  equality  constraints  (as  used  for  the  end  condition  for  the  maneuver).  For  the 
supercavitating  vehicle  system,  the  penalty  multipliers,  Pme?  and  Peq  are  positive, 
and  these  multipliers  multiply  all  equality  constraints  by  the  same  factor,  and  all 
inequality  constraints  by  the  same  factor;  in  the  general  case,  these  may  be  non¬ 
linear  and  can  be  independent  for  each  constraint. 

F(ui,u2,  ...,us,T)=T  +  P ineq  Y  M  +  peqY  \ce(h\  (5-n) 

Ci>  0  j 
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Coupled  to  the  penalty  method  is  an  unconstrained  optimizer.  Simple  search 
algorithms  worked  best  for  these  complicated  systems,  with  the  simplest  built  in 
Matlab  function  patternsearch  providing  the  best  results  in  the  fastest  time.  This 
observation  is  not  surprising,  since  the  gradient  based  methods  can  be  expected 
to  face  difficulties  due  to  the  non-smooth  dynamics,  and  the  overall  complexity  of 
the  optimization  function.  Additionally,  these  gradients  are  not  readily  available 
due  to  the  numeric  nature  of  the  dynamics  integration  and  the  cavity  (planing 
force)  calculations.  Overall,  all  algorithms  using  the  penalty  method  provided  better 
improvement  on  initial  guesses  then  the  constrained  approach  with  fmincon. 

By  using  the  penalty  approach,  the  obstacle  scenarios  similar  to  the  cylindrical 
case  shown  in  Figure  5.17,  can  be  considered  for  the  non-cylindrical  case.  The  ma¬ 
neuver  consisted  of  an  end  condition  of  (zf,Xf)  =  (—20.0 m,  80.0m)  with  an  obstacle 
centered  at  ( z0bst,x0bst )  =  (—5.0m,  40.0m).  Medium  radius  obstacles  are  solved  ini¬ 
tially,  and  these  solutions  are  used  to  seed  optimization  runs  for  progressively  larger 
radius  obstacles  until  no  feasible  solution  could  be  found.  A  feasible  solution  for  the 
non-cylindrical  single  obstacle  case  is  shown  in  Figure  5.24  for  an  obstacle  radius 
robst  =  7.0  m.  In  this  plot,  the  vehicle  and  cavity  orientation  plots  are  exaggerated 
in  size.  When  compared  to  the  cylindrical  case,  the  non-cylindrical  case  was  also 
found  to  be  capable  of  maneuvering  around  similarly  sized  obstacles. 

The  Matlab  code  for  carrying  out  the  optimal  control  with  the  non-cylindrical 
planing  model  is  presented  in  Appendix  B.  This  is  included  to  demonstrate  the  gen¬ 
eral  method  for  how  the  optimization  using  patternsearch  with  the  penalty  method 
is  coded,  along  with  the  specific  code  for  the  non-cylindrical  planing  dynamics. 
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Figure  5.24:  Non-cylindrical  planing  model,  run  for  (zf,Xf)  =  (—20.0m,  80.0m), 
obstacle  at  ( z0bst,x0bst )  =  (— 5.0m,  40.0m),  and  r0bst  =  7.0  m. 

5.3.2  Bootstrapping  Techniques  Using  Simple  Integration 

The  use  of  simple  integration  schemes  and  bootstrapping  techniques  are  com¬ 
mon  for  numeric  optimal  control  solutions.  Simple  integration  schemes,  such  as  mid¬ 
point  integration  shown  in  Eq.  (5.12),  help  speed  up  the  integration  of  the  dynamics, 
which  in  turn  dramatically  increases  the  speed  of  the  optimization.  Bootstrapping 
is  based  on  building  solutions  from  previous  solutions.  A  common  approach  for  solv¬ 
ing  a  complicated  optimal  control  problem  is  to  first  solve  the  maneuver  by  using 
a  coarse  discretization  using  a  simple  integration  scheme.  Since  the  discretization 
is  coarse,  and  the  integration  method  is  simple,  this  initial  problem  is  not  difficult 
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to  solve.  This  solution  can  then  be  applied  to  a  finer  discretization  and  solved 
again.  This  can  be  carried  out  for  several  iterations,  and  as  the  time  discretization 
becomes  finer,  the  solution  becomes  better,  and  the  integration  becomes  more  accu¬ 
rate.  Although  the  optimizer  has  to  solve  a  series  of  progressively  more  complicated 
problems,  it  is  given  a  series  of  progressively  better  starting  guesses. 


xi+i-Xi=  f  +  f(t,  x(t))dt  »  /(-  +  *t+1 ,  —  +  'I  t+1  )(b:+ 1  -  U)  (5.12) 

Ju  2  2 

This  technique  has  been  applied  to  supercavitating  vehicle  systems  before 
[39,  1];  however,  the  solutions  for  the  maneuvers  considered  are  all  within  the  cav¬ 
ity.  For  the  maneuvering  studies  in  this  work,  and  the  parameter  values  and  speeds 
considered,  planing  or  cavity  contact  is  abundant  during  maneuvering  (operation 
in  ranges  where  unstable  and  limit-cycle  behavior  is  observed  with  the  uncontrolled 
and  feedback  controlled  systems).  Bootstrapping  schemes  and  mid-point  integration 
were  attempted  for  both  the  set  end  time  and  floating  end  time  maneuvers.  The 
system  dynamics  for  non-smooth  systems  can  very  greatly  depending  on  the  re¬ 
gion  of  operation.  The  supercavitating  system  is  an  example  of  this;  when  planing, 
the  system  dynamics  are  dramatically  different,  with  much  higher  forces  present. 
Since  planing  is  prevalent,  mid-point  integration  for  any  coarse  discretization  be¬ 
comes  highly  inaccurate.  Feasible  solutions  for  coarse  discretizations  where  planing 
is  present  become  extremely  difficult  to  resolve.  Even  when  feasible  solutions  can 
be  solved,  the  dynamics  are  so  dependent  on  the  region  of  operation  that  differ¬ 
ent  solutions  “bootstrapped”  to  finer  discretizations  are  not  necessarily  feasible,  let 
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alone  good  seeds.  Because  of  this,  for  the  supercavitating  vehicle  systems,  direct 
integration  methods  are  required  for  accuracy.  This  of  course,  comes  at  the  expense 
of  lengthy  computational  time  for  the  optimization  process. 

5.4  Maneuvering  for  Delayed  Model  without  Steady  Planing  As¬ 
sumption 

The  penalty  approach  was  also  applied  to  the  delayed  model  presented  in  the 
previous  chapter.  This  model  includes  the  forces  due  to  the  motion  of  the  body, 
instead  of  the  use  of  the  steady  planing  assumption.  Due  to  the  delay,  the  equations 
of  motion  must  be  integrated  by  using  a  delay  differential  equation  (DDE)  solver  that 
requires  the  state  history  along  with  an  initial  condition.  In  the  previous  examples, 
the  equations  of  motion  are  integrated  over  each  individual  section  of  piecewise 
constant  control.  This  becomes  difficult  to  implement  with  the  DDE  solver,  since 
each  individual  section  will  require  a  state  history. 

To  get  around  this  difficulty,  a  smooth  control  scheme  is  implemented  by  us¬ 
ing  a  splined  interpretation.  In  this  formulation,  the  optimization  variables  become 
inputs  to  a  spline  interpolation  that  is  used  as  the  control  input  function;  the  op¬ 
timization  variables  in  this  case  can  be  though  of  as  parameters  in  a  describing 
function  that  defines  the  controller  input.  Now,  the  control  input  function  for  the 
duration  of  the  maneuver  is  known  completely  within  the  integration  function.  An 
additional  advantage  in  this  approach  is  that,  the  controller  input  is  now  a  con¬ 
tinuous  and  differentiable  function  which  is  practical  to  implement.  Additionally, 
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control  contributions  to  the  system  dynamics  will  be  less  abrupt.  It  should  be  noted 
that  the  previous  systems  can  be  also  solved  in  a  similar  manner,  but  in  these  cases, 
the  piecewise-constant  control  is  sufficient  to  demonstrate  results. 

Results  for  a  move-to-point  maneuver  to  ( Zf,Xf )  =  (—20.0m,  80.0m)  and  the 
delayed  model  with  non-steady  planing  are  shown  in  Figure  5.25.  The  smooth 
outer-loop  control  inputs  are  shown  in  Figure  5.26.  Obstacle  avoidance  maneuvers 
are  then  considered  starting  with  maneuvers  with  small  obstacles  and  progressing  to 
larger  obstacles  until  no  feasible  solutions  can  be  found.  The  ending  largest  feasible 
obstacle  radius  of  r0bst  =  5.0  m  for  a  maneuver  to  ( Zf,Xf )  =  (—20.0m,  80.0m)  with 
the  obstacle  centered  at  ( z0bst ,  x0bst )  =  (—5.0m,  40.0m)  is  a  smaller  radius  obstacle  as 
compared  to  similar  runs  for  the  cylindrical  and  non-cylindrical  models  with  steady 
planing  and  without  time  delay.  The  best  solved  trajectory  is  shown  in  Figure 
5.27.  Here,  it  can  be  seen  that  the  trajectory  moves  under  the  obstacle  (in  terms  of 
positive  z),  as  opposed  to  over  as  with  the  steady  planing,  non-dclayed  cylindrical 
and  non-cylindrical  models. 

The  influence  of  the  initial  guess  on  the  ending  solution  can  be  seen  by  consid¬ 
ering  a  larger  obstacle  (infeasible  ending  trajectory).  A  run  using  the  best  discovered 
trajectory  at  rQbst  =  5.0  m  as  an  initial  guess  for  the  optimization  to  do  a  maneuver 
with  an  obstacle  radius  r0bst  =  6.0  m  is  shown  in  Figure  5.28.  Here,  the  optimization 
scheme  considers  trajectories  under  (with  respect  to  positive  z)  the  obstacle,  which 
is  the  case  for  the  trajectory  that  seeded  the  optimization.  If  instead  zero  outer-loop 
control  input  is  used  as  the  initial  guess,  the  optimization  considers  trajectories  over 
the  obstacle  (see  Figure  5.29).  Although  the  results  of  either  seeding  in  this  case 
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Figure  5.25:  Delayed  model  without  steady  planing  assumption  and  run  to  (. Zf,Xf ) 


(—20.0 m,  80.0m). 
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Figure  5.26:  Outer-loop  control  inputs  for  delayed  model  without  steady  planing 
assumption  and  run  to  ( Zf,Xf )  =  (—20.0m,  80.0m). 
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Figure  5.27:  Delayed  model  without  steady  planing  assumption,  run  for  ( Zf,Xf ) 
(—20.0m, 80.0m),  obstacle  at  (z0bst,  x0bst)  =  (—5.0m, 40.0m),  and  robst  =  5.0  m. 


112 


are  infeasible,  the  difference  in  the  types  of  trajectories  considered  demonstrate  the 


influence  of  the  initial  guess. 


Figure  5.28:  Delayed  model  without  steady  planing  assumption,  run  for  ( Zf,Xf )  = 
(—20.0m,  80.0m),  obstacle  at  (z0bst,x0bst)  =  (—5.0m,  40.0m),  and  robst  =  6.0  m.  The 
best  found  control  input  solved  for  a  smaller  obstacle  run  is  used  as  an  initial  guess. 
Infeasible  ending  trajectory. 


The  Matlab  code  for  carrying  out  the  optimal  control  with  the  delayed  model 
is  presented  in  Appendix  B.  This  is  included  to  demonstrate  the  general  method  for 
how  the  optimization  using  the  smooth  inputs  is  coded,  along  with  the  specific  code 
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Figure  5.29:  Delayed  model  without  steady  planing  assumption,  run  for  ( Zf,Xf )  = 
(—20.0m, 80. 0m),  obstacle  at  ( z0bst,x0bst )  =  (— 5.0m, 40.0m),  and  r0bst  =  6.0  m. 
Trivial  outer-loop  control  is  used  as  an  initial  guess.  Infeasible  ending  trajectory. 
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for  the  impacting  and  delayed  dynamics. 


5.5  Maneuvering  with  Six  DOF  Model 

The  optimal  control  solution  method  is  also  applied  to  a  six  degree-of-freedom 
model.  The  model  is  used  to  consider  more  general  flight  motions  in  three  dimen¬ 
sions,  and  the  details  are  presented  in  Appendix  A.  This  model  does  not  include 
small  angle  assumptions  for  the  vehicle  motions,  although  there  are  assumptions  for 
the  forces  generated  by  the  control  elements.  Initially,  only  dive-plane  maneuvers 
are  considered  to  test  the  approach.  When  considering  a  fixed  end  time  maneuver, 
the  system  was  found  to  suffer  from  several  local  minima.  By  using  a  constrained 
optimizer,  it  was  found  that  for  a  maximum  angle  dive  or  climb  maneuvers  (max¬ 
imizing  or  minimizing  Of),  the  optimization  would  not  progress  far  from  an  initial 
guesses  before  becoming  stuck  in  a  local  minimum.  More  difficult  maneuvers  such 
as  move  to  point  were  difficult  to  solve,  even  when  given  a  near  feasible  starting 
guess.  Again,  a  penalty  method  based  on  the  pattern  search  algorithm  was  applied 
instead  and  was  found  to  worked  well.  A  solution  for  a  move  to  point  maneuver 
with  ( Zf,Xf,yf )  =  (—50.0m,  120.0m,  0.0m)  is  shown  in  Figure  5.30.  Feedback  con¬ 
trol  of  the  form  shown  in  Eq.  (5.9)  is  utilized  to  stabilize  vehicle  motions  in  the 
vertical  plane,  with  kinner  =  0.7.  The  vehicle  parameters  remained  consistent  with 
the  previous  simulations  with  g  =  9.81  m/s2,  m  =  22.7005  kg,  Rn  =  0.0191  rn, 
R  =  0.0508  m,  L  =  1.8  m,  and  Cx o  =  0.82.  The  propulsive  force  is  set  as  constant, 
with  Fprop  =  2200  N.  The  initial  conditions  are  for  a  straight  and  level  flight  with 
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Figure  5.30:  Six  degree-of- freedom  model,  run  for  ( Zf,Xf,yf )  =  (—50 m,  120m,  0m). 


a  forward  velocity  of  u  =  75  m/s. 

With  the  six  degree-of- freedom  model,  motions  outside  of  the  vertical  plane 
can  also  be  considered.  To  allow  for  motion  in  the  horizontal  plane,  the  optimizer 
is  allowed  three  control  inputs,  the  cavitator  actuation  angle,  the  elevator  actuation 
angle,  and  the  rudder  actuation  angle.  Both  elevators  and  both  rudders  are  assumed 
to  move  together  with  the  same  actuation  angle.  A  move-to-point  maneuver  for 
(zf,Xf,yf)  =  (—50.0m,  120.0m,  5.0m)  is  shown  in  Figure  5.31.  The  resulting  best 
solution,  is  not  as  expected  with  an  abrupt  maneuver  in  the  horizontal  plane  only 
within  the  last  time  segment  of  the  optimization  discretization.  This  observation 
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is  confirmed  by  looking  at  the  resulting  control  input,  shown  in  Figure  5.32,  where 
the  rudder  input  is  only  applied  during  the  last  time  segment.  The  result  can  be 
traced  to  the  fact  that  the  inner-loop  controller  is  only  stabilizing  planing  motions 
in  the  vertical  plane.  This  is  in  part  due  to  the  orientation  of  the  cavitator  actuation 
angle  that  is  set  to  be  vertical  (with  respect  to  the  body  orientation).  With  this, 
the  cavitator  angle  control  actuation  of  an  un-rolled  vehicle  cannot  be  expected  to 
produce  contributions  in  the  horizontal  plane.  The  vehicle  motions  in  the  horizontal 
plane  are  in  fact  quite  unstable,  with  significant  planing  forces  being  generated 
during  the  quick  motion  at  the  end  of  the  maneuver,  which  is  required  to  achieve 
desired  horizontal  displacement.  If  the  vehicle  were  to  generate  horizontal  motions 
early  in  the  maneuver,  the  high  planing  forces  and  instabilities  would  greatly  reduce 
the  ability  of  the  vehicle  to  reach  the  desired  end  point. 

The  maneuvering  studies  for  the  six  degree-of- freedom  model  are  very  prelim¬ 
inary.  More  work  needs  to  be  done  such  as  allowing  greater  control  authority  over 
vehicle  roll,  by  giving  independent  control  over  each  set  of  rudders  and  elevators. 
Better  control  of  vehicle  roll  may  allow  for  inner-loop  feedback  stabilization  in  mul¬ 
tiple  planes  with  a  single  axis  cavitator.  Additionally,  a  two  axis  cavitator  could 
also  be  considered,  with  separate  feedback  control  along  each  direction. 

5.6  Discussion 

Maneuvering  of  non-smooth  vehicle  systems  is  considered  in  this  chapter  with 
the  primary  example  of  supercavitating  vehicle  systems.  Maneuvers  were  solved  for 
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Figure  5.31:  Six  degree-of-freedom  model,  run  for  ( Zf,Xf,yf ) 

(— 50.0m,  120.0m,  5.0m). 


the  cylindrical,  non-cylindrical,  six  degree-of-freedom,  and  the  delayed  non-steady 
planing  models  described  in  the  previous  chapters. 

In  this  work,  the  speeds  and  parameter  values  for  the  maneuvers  considered 
for  the  supercavitating  vehicle  systems  are  within  the  range  where  unstable  as  well 
as  limit-cycle  behavior  exists  in  the  un-controlled  and  feedback  systems.  Within 
these  parameter  ranges,  the  non-smoothness  of  the  planing  must  be  addressed  in 
the  optimal  control  solution  procedure,  and  planing  is  apparent  in  the  best  calcu¬ 
lated  maneuvers.  As  such,  the  need  for  an  inner-loop  feedback  controller  is  required. 
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Figure  5.32:  Outer-loop  control  inputs  for  six  degree-of-freedom  model,  run  for 
(zf,Xf,yf)  =  (—50.0m,  120.0m,  5.0m). 

This  fast  acting  stabilizing  inner-loop  controller  is  used  in  conjunction  with  a  coarser 
outer-loop  control  that  guides  the  vehicle  through  the  desired  maneuver.  Addition¬ 
ally,  since  the  vehicle  moves  across  regions  with  dramatically  different  dynamics, 
simple  integration  techniques  with  coarse  time  discretizations  do  not  accurately  es¬ 
timate  vehicle  motions.  Therefore  methods  of  bootstrapping  (in  terms  of  using 
simple  integration  techniques  and  moving  from  coarse  to  finer  time  discretizations) 
are  not  viable  for  solving  maneuvers  for  these  types  of  systems. 
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For  systems  with  complicated  dynamics,  constrained  optimizers  can  face  diffi¬ 
culties  searching  for  feasible  optimal  solutions.  Penalty  methods  open  up  the  feasi¬ 
ble  domain  and  allow  for  simultaneous  progression  towards  better  and  more  feasible 
solutions.  The  penalty  approach  combined  with  a  simple  unconstrained  search  algo¬ 
rithms  (such  as  patternsearch)  where  shown  to  provide  a  means  to  solve  difficult  ma¬ 
neuvers  for  systems  with  complicated  non-smooth  dynamics  (where  the  constrained 
optimizers  failed).  Seeding  algorithms  such  as  those  for  solving  progressively  more 
aggressive  maneuvers  (iterating  towards  larger  obstacles)  using  previous  solutions  as 
initial  guesses  can  help  to  speed  up  the  optimization  process.  However,  the  influence 
of  the  initial  guess  must  also  be  considered  when  using  these  approaches,  as  they 
may  limit  the  types  of  trajectories  considered.  It  may  be  beneficial  to  use  different 
initial  guesses  once  the  iterations  terminate  with  infeasible  solutions. 
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Chapter  6 


Summary  and  Recommendations  for  Future  Work 
6.1  Summary 

This  dissertation  work  is  centered  around  supercavitating  vehicle  systems  with 
“soft”  non-smooth  interactions  between  the  cavities  surrounding  the  vehicles  and  the 
vehicles.  Here,  the  modeled  cavity  interactions  generate  a  system  that  is  charac¬ 
terized  by  non-constant  switching  boundaries  and  non-constant  switched  dynamics. 
Since,  the  vehicle  motions  vary  greatly  depending  on  the  region  of  interaction,  this 
creates  a  very  complex  system  given  that  the  boundaries,  and  the  forces  within  each 
region,  are  state  and  control  dependent. 

Much  of  the  above  mentioned  complexity  is  due  to  the  inclusion  of  realistic 
physical  effects  such  as  the  planing  associated  with  shifted  non-cylindrical  cavities, 
an  aspect  unique  to  this  dissertation  work.  By  using  these  shifted  non-cylindrical 
cavities,  similar  qualitative  changes  were  found  as  in  previous  supercavitating  ve¬ 
hicle  studies.  Similar  stabilization  techniques  were  also  successful  for  inclusion  in 
the  systems  considered  here.  However,  the  values  of  equilibrium  solutions  and  the 
parameter  values  at  which  qualitative  changes  experienced  by  them  occur,  differ  for 
the  non-cylindrical  cavity  case.  Hence,  the  inherent  vehicle  motions  when  including 
these  effects  can  be  distinctly  different. 

Additional  observations  have  been  made  by  generating  a  preliminary  vehicle 
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motion  model  for  partially  cavitating  vehicles.  This  model  is  meant  to  be  the  first 
step  towards  modeling  full  vehicle  missions  that  include  cavity  growth  and  collapse. 
For  small  cavities,  it  has  been  found  that  the  planing  is  insufficient  to  support  the 
rear  of  the  vehicle.  Furthermore,  the  use  of  passive  fin  input  which  supports  the 
vehicle  for  the  supercavitating  case,  is  not  necessarily  capable  of  rejecting  transient 
motions  that  move  the  vehicle  away  from  a  straight  and  level  flight.  Linear  fin 
feedback  is  shown  to  work  well  in  conjunction  with  the  linear  cavitator  feedback, 
and  active  fin  inputs  may  be  necessary  for  sustaining  partial  cavity  flight. 

A  full  representation  of  the  derivation  of  the  immersion  depth,  and  immersion 
rate  terms  used  in  the  planing  force  calculation  have  also  been  given  in  this  work. 
These  terms  have  been  presented  in  a  somewhat  arbitrary  manner  in  much  of  the 
previous  research.  This  full  derivation  provides  a  basis  for  properly  accounting  for 
the  effects  of  vehicle  motions  and  body  velocities  into  the  cavity,  creating  damping¬ 
like  contributions  in  the  planing  force  formulation.  This  is  a  departure  from  the 
previous  steady  planing  assumption  based  studies,  which  form  a  vast  majority  of 
the  previous  literature;  in  them,  one  only  considers  planing  forces  due  to  the  relative 
vehicle-cavity  positions.  The  complete  representation  of  the  immersion  terms  also 
allows  for  a  proper  handling  of  the  cavity  time-delay  effect  (in  this  case  for  cylindrical 
cavities).  In  this  work,  the  delay  is  found  to  have  a  stabilizing  effect  for  particular 
values  of  cavitation  numbers. 

A  combined  inner-loop  and  outer-loop  control  scheme  is  applied  with  suc¬ 
cess  for  the  maneuvering  studies.  Here,  fast  acting  instabilities  are  rejected  by  a 
feedback  inner-loop  while  a  numeric  optimal  control  derived  outer-loop  guides  the 
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vehicle  through  a  desired  maneuver.  For  complicated  supercavitating  vehicle  system 
maneuvers,  the  best  results  are  provided  by  utilizing  a  penalty  method  (to  account 
for  the  constraints)  along  with  a  simple  unconstrained  search  algorithm.  This  is  due 
to  the  highly  complicated  feasible  domain  generated  by  dynamic  constraints.  Direct 
integration  techniques,  rather  then  boot-strapping  techniques  are  required,  since  the 
system  dynamics  differed  greatly  depending  on  operating  conditions.  Maneuvers  are 
generated  for  cylindrical,  shifted  non-cylindrical,  six  degree-of-freedom  systems,  and 
impacting  models.  All  maneuvers  have  been  performed  at  speeds  where  there  were 
tight  cavity-body  clearances,  and  planing  is  dominant  during  the  motion.  Much  of 
this  work  can  be  extended  for  use  with  other  non-smooth  systems. 

6.2  Recommendations  for  Future  Work 

There  are  many  direct  extensions  possible  from  this  body  of  work.  Full  mis¬ 
sion  simulation  by  using  the  partial  cavitation  model  is  still  an  open  avenue.  The 
main  limitation  for  this  line  of  work  is  the  availability  of  accurate  partial  cavitation 
models,  which  in  particular  should  include  non-steady  cavity  growth  and  collapse. 
The  numeric  optimal  control  approach  outlined  for  maneuvering  should  work  well 
for  developing  optimal  outer-loop  commands  for  these  full  mission  simulations.  The 
six  degree-of-freedom  model  is  also  another  direction  for  expansion.  Out  of  plane 
maneuvers  with  multiple  axis  cavitator  actuation  can  be  considered.  The  non- 
cylindrical  cavities  accounting  for  delay  affects  have  also  not  been  considered.  This 
may  greatly  complicate  the  computation,  since  multiple  delays  accounting  for  cav- 
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ity  positions  all  along  the  planing  area  must  be  tracked.  Additional  physical  effects 
such  as  gravity  effects  on  the  cavity  should  also  be  considered,  since  as  shown,  the 
cavity  shape  can  play  a  significant  role  in  determining  the  system  dynamics.  The 
approach  used  to  determine  the  non-cylindrical  planing  forces  in  this  work  is  gen- 
eralizable  to  other  cavity  models  that  are  to  account  for  realistic  physical  effects. 
Actuator  dynamics,  specifically  actuator  rate  limitations,  has  not  been  included  in 
this  work,  but  can  be  incorporated  into  the  numeric  optimal  control  approach.  The 
maneuvering  studies  can  also  be  extended  to  a  broad  range  of  non-smooth  vehicle 
systems,  including  hypersonic  flight  vehicles. 
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Appendix  A 


Six  Degree-of-Freedom  Model 

A  six  degree-of- freedom  (DOF)  model  is  generated  to  consider  general  flight 
dynamics  modeling.  In  this  model,  all  small  angle  assumptions  related  to  the  dy¬ 
namics  are  removed,  and  the  vehicle  speed  is  no  longer  considered  to  be  constant. 
For  supercavitating  vehicle  systems,  as  with  many  flight  systems,  the  force  defini¬ 
tions,  position  vectors,  and  control  surface  rotations  are  conveniently  described  in 
local  coordinate  systems.  An  Euler  angle  approach  is  utilized  to  track  several  coor¬ 
dinate  systems  and  their  relationships  to  an  inertial  reference  frame.  This  approach 
closely  mimics  the  one  used  in  previous  literature  [13]. 

A.l  General  Approach 

An  inertial  reference  system  with  made  up  of  unit  vectors  <  ehe^e^  >  is 
defined.  A  moving  reference  frame  that  is  attached  to  the  body  CG  and  aligned 
along  the  body  is  defined  with  unit  vectors  <  bi,  b-2 ,  63  >,  as  shown  in  Figure  A.l.  By 
using  the  body  coordinate  system,  the  locations  to  the  control  surfaces  are  simply 
defined.  Additionally,  the  vectors  defining  the  forces  and  moments  can  easily  be 
resolved  and  applied  within  this  reference  system  (similar  to  the  moving  coordinate 
system  in  the  dive-plane  model). 

The  Euler  angle  relationships  between  the  differing  coordinate  systems  used 
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Figure  A.l:  Diagram  of  inertial  and  body  reference  frames. 

in  this  work  follow  the  ZYX  convention.  An  example  of  transforming  from  the 
inertial  frame  <  ei,e2,e3  >,  to  the  body  frame  <  bi,b2,b3  >,  is  shown  in  Figure 
A. 2.  The  first  rotation  is  about  63  axis  with  angle  T,  generating  the  intermediate 
frame  <  i\  ,e2  ,63  >.  The  second  rotation  is  about  b2  with  angle  0,  generating  the 
intermediate  frame  <  e\",  e2" ,  £3  >.  The  final  rotation  is  about  e  ”  with  angle 
to  the  body  fixed  frame.  The  rotation  matrix  [Re->b\  that  dehnes  the  transformation 
from  the  inertial  to  body  reference  frame  can  then  be  expressed  as  a  chain  of  rotation 
matrixes  as  shown  in  Eq.  (A.l). 


126 


<I> 


Figure  A. 2:  ZYX  rotation  order  used  for  Euler  angle  relationships. 
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The  cavitator  reference  frame  is  one  that  follows  the  orientation  of  the  cavitator 
as  shown  in  Figure  A. 3.  The  cavitator  is  constrained  to  only  one  control  actuation 
angle  5C  about  the  b2  or  c2,  axis.  So  the  cavitator  reference  frame  can  be  described 
by  a  single  rotation  from  the  body  reference  frame  as  shown  in  Eq.  (A. 2). 
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The  fin  reference  frames  can  be  represented  as  shown  in  Figure  A. 4.  For  each 
fin,  f\  is  oriented  forward  along  the  fin,  and  /2  extends  out  along  the  hn  axis.  Two 
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Figure  A. 3:  Diagram  of  cavitator  reference  frame. 
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Figure  A. 4:  Diagram  of  fin  reference  frames. 

successive  rotations  define  the  fin  reference  frame,  first,  a  sweepback  angle  rotation 
about  the  positive  /3  axis,  and  a  control  angle  5f  rotation  along  the  fin  axis,  /2.  Each 
fin  is  associated  with  a  unique  transformation  between  the  body  and  fin  reference 
frames. 
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Figure  A. 5:  Velocity  direction  shown  with  respect  to  a  local  coordinate  system 

<  Xi,X2,X3  >. 

A. 2  Fin  Forces 

The  velocity  at  any  point  can  be  calculated  by  using  rigid  body  kinematics. 
This  velocity  can  then  be  expressed  in  any  of  the  local  reference  frame  basis  by 
using  the  proper  transformation  matrixes.  An  example  of  a  velocity  vector  shown 
with  respect  to  some  generic  reference  system,  <  Xi,x2,x3  >,  is  shown  in  Figure 
A. 5.  The  velocity  can  be  resolved  in  the  local  reference  frame  unit  vecotrs  as  V  = 
ux\  +vx\  +vjx\  .  The  individual  slip  angles  a  and  /3  can  then  be  expressed  according 
to  Eqs.  (A.3)-(A.4).  Furthermore,  a  transformation  from  the  local  reference  frame 
to  a  frame  aligned  along  the  velocity  direction,  can  be  described  by  two  rotations 
about  the  slip  angles  a  and  /3.  The  transformation  matrix  can  be  described  by  Eq. 
(A. 5).  The  slip  angles  are  important  in  determining  the  forces  generated  by  the 
control  surfaces. 
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Figure  A. 6:  Velocity  direction  at  a  fin. 
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The  fin  force  is  determined  by  using  a  simplified  calculation.  The  fins  are 
modeled  as  only  generating  lift  and  drag  forces.  A  fin  diagram  and  it’s  local  velocity 
vector  are  shown  in  Figure  A. 6.  The  velocity  at  the  fin  can  be  calculated  in  the 
body  reference  frame  according  to  Eq.  (A. 6),  as  the  sum  of  the  velocity  due  to  the 
motion  of  the  CG,  along  with  the  velocity  due  to  rotation  of  the  body  (Pj  represents 
the  vector  to  the  fin  from  the  body  CG  in  the  body  reference  system).  The  velocity 
in  the  fin  reference  system  can  then  be  represented  by  a  transformation  as  shown  in 
Eq.  (A.  7). 
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(A. 6) 


v,b  =  v, l  +  ab  x  f* 

vj  =  (A. 7) 

From  the  expression  for  Vj  ,  the  slip  angles  a  and  /3  can  be  calculated  from 
Eqs.  (A.3)-(A.4),  and  the  rotation  matrix  [Rg^f]  can  be  generated.  The  lift  and 
drag  forces  are  taken  as  only  functions  of  the  angle  of  attack  a  and  they  are  defined 
according  to  Eqs.  (A.8)-(A.9).  The  lift  and  drag  forces  for  a  control  surface  are 
generally  defined  with  respect  to  the  relative  fluid  velocity  direction  as  shown  in 
Figure  A.  7.  Therefore  the  lift  and  drag  forces  can  be  directly  represented  in  the 
velocity  reference  frame,  and  the  overall  force  of  the  fin  can  be  expressed  according 
to  Eq.  (A.  10),  when  transformed  back  into  the  body  reference  frame. 
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A. 3  Cavitator  Force 


(A.10) 


The  cavitator  force  is  calculated  somewhat  differently  due  to  the  cavitation.  A 
diagram  of  the  orientation  of  the  lift  and  drag  forces  for  the  cavitator  with  respect  to 
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Figure  A.  7:  Lift  and  drag  forces  on  a  control  surface  with  respect  to  the  relative 
fluid  velocity  direction. 

the  relative  fluid  velocity  direction  are  shown  in  Figure  A. 8.  Through  experimental 
data,  for  a  cavitating  disc,  it  was  found  that  the  forces  on  the  cavitator  were  pre¬ 
dominantly  due  to  pressure  [18].  This  means  that  the  lift  to  drag  ratio  follows  the 
relationship  fij fd  ~  tan(ac),  or  that  the  force  is  predominantly  along  the  direction 
perpendicular  to  the  wetted  face  of  the  cavitator.  From  experimental  results,  the 
coefficient  of  drag  for  flow  along  the  axial  direction  of  a  disc  cavitator  is  found  to 
follow  a  relationship  according  to  Eq.  (A.  11).  So  the  overall  cavitator  force  can  be 
represented  as  given  in  Eq.  (A.  12),  where  cos(ac)  =  pp.  The  cavitator  force  fc  is 
applied  along  the  axial  direction  of  the  cavitator,  so  the  vector  representation  of  the 
cavitator  force  in  the  body  reference  system  can  is  expressed  as  in  Eq.  (A.  13). 
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Figure  A. 8:  Lift  and  drag  forces  for  a  cavitator  with  respect  to  the  relative  fluid 
velocity  direction. 
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A. 4  Cavity  and  Planing  Force 


Since  various  speeds  are  considered  in  the  six  DOF  model,  a  closed  form  cav¬ 
ity  model  is  utilized  for  faster  generation  of  cavity  shapes.  A  cylindrical  cavity 
assumption  (without  cavity  shift  effects)  is  used  for  simplicity.  The  cavity  shape 
is  predicted  by  using  a  semi-empirical  closed  form  solution  formulated  in  reference 
[31].  As  presented  in  Chapter  3,  the  cavity  radius  at  a  point  along  the  cavity  can 
be  estimated  by  using  Eq.  (A.  14).  Here,  again  x  represents  the  length  from  the 
cavitator,  and  dc  represents  the  cavitator  diameter.  An  entire  cavity  profile  can  be 
generated  by  evaluating  Rc  at  several  points.  Alternatively,  as  carried  out  here,  an 
approximate  cylindrical  cavity  can  be  generated  by  utilizing  the  cavity  radius  at 
x  =  L. 
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dmax  =  dc\J  0.82(1  +  a)/a 
lm  =  dc/2{1.92/a  -  3) 
ki  =  1.92(0.82(1  +  a)/a)"5 
/^2  *  dc  dc)  j lm 

Rc(x)  =  dmax/2^1-(l-kl)\l-k2\y^  (A.  14) 

The  planing  formulations  presented  in  the  previous  chapters  are  for  planing  in 
two  dimensions.  Since  the  body  and  cavity  are  both  bodies  of  revolution,  and  both 
axes  share  a  common  point  (the  nose  of  the  vehicle),  the  planing  force  can  still  be 
considered  in  two  dimensions  if  an  appropriate  plane  of  planing  is  chosen.  The  plane 
of  planing  is  the  plane  in  which  the  vehicle  immersion  (when  present)  is  symmetric, 
and  hence,  the  two  dimensional  planing  models  are  valid.  This  is  the  plane  defined 
by  the  cavity  and  body  axis  (or  the  velocity  direction  and  body  axis).  A  diagram 
illustrative  of  this  relationship  is  shown  in  Figure  A. 9.  Here,  the  rotated  planing 
frame  <  Pi,p2,P3  >  is  defined  with,  pi  =  b\,  p2  =  ,  and  p3  =  pi  x  p2.  The 

planing  force  is  along  the  p3  direction  with  Fpiane  =  fpp3,  where  fp  is  the  planing 
force  calculated  by  using  the  cylinder-on-cylinder  Paryshev  formulation  presented 
in  Eq.  (2.21). 
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Figure  A. 9:  Plane  of  planing  defined  by  cavity  axis  (or  velocity  direction)  and  body 
axis. 


A.  5  Equations  of  Motion 

The  system  has  twelve  states.  The  states  considered  are  the  body  rotation 
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local  angular  velocity  is  ujb,  and  can  be  expressed  in  terms  of  the  body  rotation  angle 
rates  as  given  in  Eq.  (A.  16).  By  utilizing  transformations  from  the  intermediate 
reference  frames,  the  relationship  between  the  body  rotation  angle  rates  and  the 
angular  velocities  can  be  generated  as  shown  in  Eq.  (A.  17). 
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(A.  17) 


By  applying  Newton’s  and  Euler’s  principles,  the  sum  of  the  forces  and  the 


sum  of  the  moments  can  be  expressed  as  shown  in  Eqs.  (A.  18)- (A.  19). 
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(A. 19) 


The  forces  and  moments  that  act  upon  the  body  come  from  the  planing,  the 


four  fins,  the  cavitator,  gravity,  and  the  propulsive  force  Fprop. 
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Appendix  B 


Matlab  Code 

B.l  Partially  Cavitating  Vehicle  Dynamics 

%2/9/2009 

%This  code  is  setup  so  that  the  matrix  algebra  for  the 
equations  of  motion 

%are  done  within  the  odefunction.  This  is  required  for  the 
partial  cavity 

%simulations  since  the  added  mass  terms  cannot  be  simply 
added  as  forces 

%since  they  affect  acceleration  (inertia  matrix  becomes 
dynamic )  . 

clear  ; 

close  all  ; 

clc  ; 

global  R  L  Rn  LI  V  M  xcg  Iyy  rho  T  n  g  gcounter 
individualshift  ycptrunc  xcptrunc  Lc  rl  ; 

gcounter  =0; 
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Wo%%%%%%mODY/RUN  PARAMETERS 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

g  =  9.81; 

m=  2;  %  density  ratio 

Rn  =  0.0191; 

R  =  0.0508; 

L  =  1.8; 

sigma  =  0. 066925; 

%sigma  =  0 . 046 ;  %or  0.043  for  using  scax  (supercav)  t 
partial  cavity  shape 

11=  0.5; 

CxO  =  0.82; 

Cx  =  CxO*(l  +  sigma)  ; 

V  =  sqrt  (0.03*75~2/sigma)  ; 

rho  =  1000; 

Ll=L/3;  %olength  of  the  conical  section  of  the  body 
L2=L— LI ;  %olength  of  cylindrical  section  of  body 


create 


138 


rl=0;  %radius  of  the  front  of  the  body 


M=m*rho*pi  *  (R~2*L2+(IT3—  r  1  "3)  /  (3*  (R— r  1 ) )  *L1 )  ; 
xcg  =  1/ 4*(  — 3*R~  2* LI ~ 2+2*R*Ll ~  2*  r  1+r  1  ~  2* LI ~  2+6*R7  2*L  2 )  /  ( 3*R 
~ 2*L— 2*R  2* Ll+Ll *R* r  1+L1  * r  1  "2)  ; 
xcg=-xcg ; 

Iyy  =  l/60*m*rho  *  (3*R~4+3*r  1  *R73+12*R7  2* LI  ~2+3*R7  2*  r  1  ~2+6*R*r  1 
*L1 ~2+3*R*r 1  ~3+2*r  1  ~2*Ll~2+3*rl  "4)  *pi*Ll  +  l/12*pi*R7  2*  (L— 
LI )  *rho*m*  (3*R~2+(L— LI )  "  2)+pi*R7  2*  (L— LI )  *rho*m*  ( 1/2  *L 
+  1/2*L1)  "  2 ; 

T=.5*rho*pi*Rn~2*V*Cx; 


Wo%%%%%%%%%%%°/(CA  VITY  SHAPE  MODEL 


%data  for  initial  cavity  length  vs.  cavitation  number  to  be 
converged  upon 
%for  scax  model 

sigmafit  =[0.03649  0.039992  0.042468  0.034994  0.03196 
0.029983  0.027951  0.025967  0.024977  0.045]; 
lengthfit  =[50.934422  45.790226  42.768097  53.378616  59.052101 
63.461826  68.921989  74.727722  78.073776  40.031219]; 
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cavlen=spline  (  sigmafit  ,  lengthfit  ,  sigma); 


ficl=fopen(  ’body .  dat  ‘  ,  'w ' )  ; 
fprintf(fid  ,  ’ body  .  dat \nCONE\nC\nD\n  ’ )  ; 
fprintf(fid  ,  ’32\tl6\t  — >MBOD,  JY1END\  n  ’ )  ; 

fprintf  (  fid  ,  ’  90 . 0 \  tO  .  5\ t->.HALF.CONE.ANGLE.  (DEGREES)  ,  .HEIGHT 
OF,  ENDPLATE\n  ’ )  ; 

fprintf(fid  ,  [  num2str  (  cavlen  )  ’ \t0 . 05\ t—>. CAVITY  .LENGTH,  . 
CAVITY .NODE.FACTOR.  (CJPT)  .TO.DEFINE.#.OF.NODES.IN.CAVITY\ 
n’]); 

fprintf  (fid  ,  ’  l\t->. ITERATE  .ON.CAVITATION  JNUMBER:  ,.1=ITERATE,  . 
0=NOT.ITERATE\  n  ’ )  ; 

fprintf(fid  ,[  num2str  ( sigma )  ’  \  t5 . 0E— 5\t— >.CAVITATION.NLOVIBER. 

TO  .BE .CONVERGED. ON,  CONVERGENCE.CRITERION\  n  ’  ] )  ; 
fprintf  (fid  ,  ’  1  \  1 0 . 5  \  t  — >.l— IF  .NON-DIMENSIONAL  .WITH.BASE. 

DIAMETR,  .BASE  .RADIUS  .OTHERWISE\  n  ’ )  ; 
fprintf(fid  ,  ’2 . 0E— 4\t— >.EPS.  (EPSILON)  .for. cavity.  length\n’)  ; 
fprintf  (  fid  ,  ’  0 \ t  — > .GRAVITY .EFFECT :  .1=ADD, ;  0  IX)  NOT  ADD\n  ’ )  ; 
fprintf  (fid  ,  ’  104.0 \t— >.FROUDE  JNfUMBER\n  ’)  ; 

fprintf  (fid  ,  ’  O\t->.ANGLE.OF^TTACK.EFFECT:  .1=ADD,  .0=DO.NOT. 
ADD\n ’ ) ; 

fprintf  (  fid  ,  ’  1  0\  t->,  ANGLE  OF  ATTACK.  ( TN .  DEGREES)  \ n  ’ )  ; 
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fprintf  (  fid  ,  ’  3 . 0E7\t— >. REYNOLDS JMUMBER’ )  ; 
fclose  (  fid  )  ; 

%! sc  ax .  exe 

%fid=f open  (’  control .  xls  ’);  %use  this  incase  of  using  scax 
fi d=fopen (  ’ bookl  .  csv  ’ )  ;  %data  from  previous  run  for  pscax 
sigma  =  0.  066925; 

Cd=str2num (fgetl(fid))  ; 

Mbod=str2num  ( fgetl(fid)); 

Mcav=str2num  ( fgetl(fid)); 
counter  =  1; 
while  ( 1 ) 

readstr=fgetl(fid)  ; 
if  readstr~=— 1 

tester  (counter  ,:)  =str2num(  readstr  )  ; 
counter=counter  +  1; 

else 

break 

end 

end 
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fclose  (  fid  )  ; 


xpos=tester ( : , 1 ) ; 
ypos=tester ( : ,2) ; 

xcp  =  [xpos  (Mbod+1) ;  tester  (Mbod+l:Mcav+Mbod, 3)  ]  ;  %adds  first 
x,y  value  to  the  cavity  section  of  xcp/ycp  data 
ycp  =  [ypos  (Mbod+1) ;  tester  (Mbod+l:Mcav+Mbod, 4)  ]; 

%xcp  =[x  (Mbod+1) ;  xcp  ] ; 

%ycp  =[y  (Mbod+1) ;  yep]; 

normLF=L/(Rn*2)  ;  %finds  the  normalized  value 

for  L 

Lc^nax(  xcp )  *R,n*  2 ; 

if  L<Lc  %oonly  truncates  if  cavity  length 

>L 

index=find  (xcp<normL ,  1,  ’last  ’  );  %finds  the  last  x 

point  index  that  is  <  L 

xcptrunc  =  [xcp  ( 1 :  index  )  ;  normL  ]  ;  %integrates  only 

up  to  length  =  L 

ycptrunc  =  [ycp  ( 1 :  index  )  ;  interpl  (xcp  ,  yep  ,  normL)  ]  ; 

else 

xcptrunc=xcp  ; 
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ycptrunc=ycp  ; 


end 


Wo%%%%%CODE  FOR  MANUALLY  MODIFIED  CAVITY  SHAPES  %necessary 
when  using 

Wo%%%%%scax  model  to  create  partial  cavity 

if  y cp tr unc (end) *Rn*2  <  R 

index=find  ( (  ycptrunc*Rn*2)>R,  1  ,  ’last  ’ )  ; 
xcptr unc=xcptrunc ( 1 : index )  ; 
yeptr unc=ycptrunc ( 1 : index )  ; 

end 

Wo%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


individualshift=— 1 Cd/8*  cumtrapz  (xcptrunc  ,  1 .  /  (  ycptrunc  . "  2)  )  ; 


Wo%%%%DDE%%%%%mTEGRA  TION 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

delta_t=0.0001;  x_init=[0.1  0.2  0  0];  i  =  1 ; 
x_init=[0  0  0  0]; 
for  tl  =0:  delta_t  :  1 

[t_state  ,  x.state]  =  ode45  ( @ODEfunPC_mxinside_020809  ,  [ 
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tl  tl+delta_t],  x_init); 


x  ( i  ,  :  )  =x_st  at  e  (end  ,:)  ;  t(l,i)=t_state  (end)  ; 
x_init=x_state  ( end  ,  :  )  ; 

%size  ( t-state  ) 

clear  x_state  t_state; 

i=i+l 

end 


% 

figure  ( 1 )  ; 

subplot  (2,2,1)  ;  plot  ( t  ,x(:  ,1)  )  ;ylabel(  ’  depth  uz  u  (m)  ’ )  ;  grid  on  ; 
xlabel  (  ’  Time^  (  s  )  ’ )  ; 

subplot  (2,2,2)  ;plot(t  ,x(:  ,2)  )  ;ylabel(  ’velocity  (m/  s  )  ’ )  ; 
grid  on  ;  xlabel  (’  Time^  (  s  )’)  ; 

subplot  (2,2,3)  ;plot(t  ,  x  ( :  ,3)  )  ;ylabel(  ’pitch  ^angle^\theta^( 
rad )  ’ )  ;  grid  on  ;  xlabel  (  ’  Time^  (  s  )  ’ )  ; 
subplot  (2,2,4)  ;  plot  ( t  ,x(:  ,4)  )  ;ylabel(  ’pitcher  a  te,^q^(rad/s)  ’) 
;  grid  on  ;  xlabel  (’ Time^  (  s  )’)  ; 
function  dxdt  =  ODEfunPC_mxinside_020809  ( t  ,  x) 

%  state  model  from  Dzielski  &  Kurdila  ’s  paper 
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%  x—[z  w  theta  q] 


global  R  L  Rn  LI  V  M  xcg  Iyy  rho  T  n  g  A  B  C  D  gcounter 
individualshift  ycptrunc  xcptrunc  Lc  rl  ; 

%%%OREA  TING  THE  MX7  s  FOR  THE  EX)1 \ff%%%%%%%%%%%%%%%%%%%%%%%%%%% 

all=T*(— 1— n)  ; 

al2=M*V— r T*n*L ; 

a21=— T*n*L ; 

a22=-M*xcg*V— ' T*n*L ~ 2 ; 


bll=— T*V*n ; 
bl2=— T*V ; 
b21=— T*V*n*L  ; 
b22  =  0; 


if  Lc<=Ll  %cavity  ends  on  forebody  (which  is 

truncated  cone) 

cl=M*g— rho*g*pi  (R"2*(L— LI)  +1/3*(L1— Lc)  *(R"3  —  ((R—  rl )  /LI* 
Lc+rl ) *  3) / (R—  ((R— rl ) /Ll*Lc+rl ) ) ) ; 
c2=-M*g*xcg— rho*g* (1 /4* pi  * (R— r 1 ) "2/Ll  *  2* (LI "4— Lc " 4) +2/3* 
pi*r 1  * (R— r 1 ) /LI  * (LI  "3— Lc " 3)  +  l/2*pi*r 1 "2*(L1"2— Lc"2) 

+  l/2*pi*R"2*(L"2— LI  "2)  )  ;  %xcg  is  negative 
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Iamda22=— (cl-M*g)  / g  ; 
lamda26=— (c2-[M*g*xcg )  /g  ; 

lamda66=rho  *  ( 1/5*  pi  *  (R— r  1 )  "  2/Ll ~  2*  (LI  "5— Lc  ~5)+l/2*pi*rl 
*(R— rl )  /  LI*  (LI  "4— Lc  ~4)+l/3*pi*rl~2*(Ll~3— Lc~3)+l/3*pi 
*R  2*(L  3 — LI  3) )  ; 

elseif  Lc<=L  %cavity  ends  on  rear  body  (cylinder) 

cl=M*g— rho*g*pi*R~2*(L— Lc)  ; 
c2=-M*g*xcg  —  ,5*rho*g*pi*R"  2*  (L"2— Lc  *  2)  ; 
lamda22=— (cl-M*g)  / g ; 
lamda26=— (c2-fM*g*xcg )  /g  ; 
lamda66  =  l/3*  pi*rho*R4  2*  (L~3— Lc  "  3)  ; 

else  %cavity  ends  past  body 

cl=M*g ; 
c2=-M*g*xcg ; 
lamda22  =  0; 
lamda26  =  0; 
lamda66  =  0; 

end 

dl  =  1; 

d2=L ; 
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%Llp=rho*  pi*R  ~2; 


MO=[M  -Al*xcg;  — M*xcg  Iyy  ]  ; 
AO=[all  al2 ;  a21  a22 ] ; 
B0=  [bll  bl2  ;  b21  b22  ]  ; 
CO=[cl  ;c2  ]  ; 

D0=[dl;d2] ; 


%Added  mass  change 


M0=  [Mflamda22  — M*xcg+lamda26 M*xcg+lamda26  Iyy+lamda66 


'o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


A2=M0\A0 ; 

B2=M0\B0 ; 

C2=M0\C0 ; 

D2=M0\D0 ; 

a22=A2  (1,1);  a24=A2  (1,2);  a42=A2(2,l);  a44=A2(2,2); 
b21=B2 (1,1);  b22=B2 (1,2);  b41=B2(2,l);  b42=B2(2,2); 
c2=C2(l,l); 

d2=D2 (1,1);  d4=D2(2,l); 
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A=[0  1  -V  0;  0  a22  0  a24 ;  0  0  0  1;  0  a42  0  a44  ]  ; 
B=  [0  0;  b21  b22 ;  0  0;  b41  b42  ]  ; 

C=  [0 ;  c2  ;  0  ;  0  ]  ; 

1)  [0 :  <12  :  0  :  <11  ] : 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


C=  [  0 ;  C2  (1)  ;0;C2(2)  ]  ;  %now  the  second  term  in  C2  is  not  =  0, 
especially  with  the  bouyancy  term,  so  it  needs  to  be 
included  unlike  before 

w  =  x(2  ; 

delta_e  =  .  12  +  .3*x  ( 4  , : )  ;  %fin  angle  wrt  body 
%delta_e=.l;%passive  fin  input 

%delta_c  =  —  15*x  ( 1  , : ) +30*x  ( 3  , : ) +0 .3*x  ( f  , : )  ;  (Kurdila  in  JVC) 
delta_c  =  15*x  ( 1  , : )  — 30*x  (3  , : )  — 0.3*x  (4  ,:)  ;%—2*x  (2  , : ) /V;  % 
correction  by  Guojian  Lin 

alphacav^w/V+delt a_c  ; 
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a  1  p  h  a— at  an  (w/V)  ; 


if  (w>0) 

cavit  y  cor  ds=abs(— ycptrunc+indi  vi  dual  shift  *sin(  alphacav)* 
cos  (  alphacav  ) )  *2*Rn;  %unnormalized  cavity  radii 
along  planing  location 

else 

cavity  cor  ds=abs  (  ycptrunc+ individuals  hi  ft  *sin(alphacav)* 
cos  (  alphacav  ) )  *2*Rn; 

end 

delt  a=abs  (  ycptrunc  *2*Rn)  — R;  %individual  deltas  ( 

difference  in  radius  of  cavity  vs  body)  unnormalized 
xcords=xcptrunc *2*Rn ;  %unnormalized  x  p  ositions 

body cor ds=abs (xcords*t an (alpha) )+R; 

%r efining  mesh  only  where  planing  begins 

index=find  ((bodycords>cavitycords)  .  *  (  xcords>max(  xcords)  /  3) 
,1);  %finds  first  location  of  planing  of  the  cylinder 
part  of  the  body 
if  index>0  %if  planing 

xcordsnew=linspace  (  xcords  ( index  —2)  ,  xcords  (end)  )  ’ ; 
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cavitycords=interpl  (xcords  ,  cavitycords  ,xcordsnew)  ; 
dclt  a=abs  ( interpl  ( xcords  ,  ycptrunc  ,  xcordsnew )  *2*Rn— R)  ; 
xcords=xcordsnew  ;  %xcords  becomes 

truncated  to  just  the  planing  area 
body cor ds=abs ( xcords  * tan ( alpha ) )+R; 

end 

hdepths  =  (bodycords— cavitycords  )  .  *  (  bodycords>cavitycords  )  ; 

%finds  immersion  depths  ( only  where  body  is  planing 

) 

if  index>0  %only  does  this  calculation  if  planing 

for  j  =1:  length  ( xcords  )  %requires  ’’for”  loop  since 

delta+hdepths  can  =0  giving  NaN 
if  ( bodycords  ( j  )>cavitycords  ( j  ) )  %calculation 
only  where  planing  (on  now  truncated  xcords) 
integrand  ( j  )  =2*  dclt  a  (j)~2./(delta(j  )+h  depths  ( j  ) ) 
3 ; 

else 

integrand  ( j  )  =0; 

end 

end 

else 
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integrand  (1:  length  (xcords)  )=0; 


end 


hdots=V*sin  ( abs  (alpha )— atan  (diff(cavitycords)  ./  diff(xcords)) 
)  ;%cavity  angle  =  atan  (  dif f  (  cavity  cor  ds  )./ diff  (  xcords  )  ) 

FpnormNoncyl=— sum(  diff  ( cumtrapz  ( xcords  ,  integrand  ’)  )  .  *  hdots 
. '  2 )  *sign  (  alpha  )  ; 

%this  DOES  take  into  account  the  cavity  slope,  d  if  f  ( cumtrapz 
)  gives  the 

%individual  trapezoidal  areas,  this  multiplied  by  hdot"2 
over  that 

%particular  area  gives  the  individual  integral  totals, 
summed  up  gives  the 

%entire  planing  force. 

if  FpnormNoncyl~=0 

xp=— sum(  diff  (  cumtrapz  ( xcords  ,  integrand  ’ .  *  xcords  )  )  .  *  hdots 
. '  2*  sign  (  alpha )  )  .  /  FpnormNoncyl ; 

else 

xp=L; 
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end 


gcounter=gcounter  +  1; 

F p=F p nor mNoncy  1  *  p  i  *  r  h  o  *R '  2 ; 

%in  original  formulation  these  terms  were  accounted  for 
within  the  D  matrix 

%and  were  therefore  omitted  from  the  planing  force  ,  in  this 
formulation  the 

%D  matrix  is  setup  to  handle  the  un—  normalized  planing  force 


dxdt  =  A*x+B*  [  delta_e  ;  delt  a_c  ]+CdD*Fp  .  *  [1  1  1  ( xp )  /L  ]  ’ ; 

B.2  Maneuvering  with  the  Non- Cylindrical  Planing  Model 

clear 

clc 

close  all 

sigma  =  0 . 03 ; 

%preamble  used  to  setup  all  parameters  needed  for  the 
function  call  .  it 
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%requires  that  sigma  be  specified  in  the  workspace  before 
calling  . 

global  Rcdot  Rp  R  V  L  A  B  C  D  Rn  Re  gcounter  individualshift 
ycptrunc  xcptrunc ; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


g  =  9.81; 
in  —  2 ; 

Rn  =  0.0191; 
R  =  0.0508; 

L  =  1.8; 


%sigma  =  0.0335;  %0.035  or  0.03 


%sigma  =0 . 025; 


%R=R* 1.25; 


%L=L  *1.1; 
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n  =  0.5; 


CxO  =  0.82; 

Cx  =  CxO*(l  +  sigma)  ; 

C  =  l/2*Cx*Rn~2/R72; 

V  =  sqrt  (0.03*75"2/sigma)  ; 

M0=  [7/9  17*L / 3 6 ;  17*L/36  1 1/60*R7 2  +  133/405*L  *  2 ] ; 

A0=C*V*[(1  —  n)  /m/L  — n/m;  —  n/m  —  n*L/m]+V*  [0  7/9;  0  17*L/36]; 
B0=C*V~2*[  —  n/m/L  1/m/L;  —n/m  0]; 

C0=  [7/9 ;  17*L / 36]  *g  ; 

D0=  [1/m/L;  1/m]; 

D00  =  [l/m/L  0;  0  1/m/L]; 

%  correction  by  Guojian  Lin 

A0=C*V*[(  —  1  —  n) /m/L  —n/m;  —  n/m  —  n*L/m]+V*  [0  7/9;  0  17*L/36] 
B0=C*V"2*[  — n/m/L  —1/m/L;  —n/m  0]; 

A2=inv  (MO)  *A0 ; 

B2=inv  (MO)  *B0 ; 

C2=inv  (MO)  *C0 ; 

D2=inv  (MO)  *D0 ; 


154 


D22=inv  (MO)  *D00 ; 


%data  for  initial  cavity  length  vs.  cavitation  number  to  be 
converged  upon 
%for  scax  model 

sigmafit  =[0.03649  0.039992  0.042468  0.034994  0.03196 
0.029983  0.027951  0.025967  0.024977  0.045]; 
lengthfit  =[50.934422  45.790226  42.768097  53.378616  59.052101 
63.461826  68.921989  74.727722  78.073776  40.031219]; 


V/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/OA 


a22=A2  (1,1);  a24=A2(l,2);  a42=A2(2,l);  a44=A2(2,2); 
b21=B2  (1,1);  b22=B2  (1,2);  b41=B2(2  ,1)  .;  b42=B2(2,2); 
c2=C2(l,l); 

d2=D2 (1,1);  d4=D2(2,l); 

d21=D22  (1,1);  d22=D22(l  ,2)  ;  d41=D22  (2  , 1 )  ;  d42=D22  (2  ,2)  ; 

A=[0  1  -V  0;  0  a22  0  a24 ;  0  0  0  1;  0  a42  0  a44  ]  ; 

B=  [0  0;  b21  b22 ;  0  0;  b41  b42  ]  ; 

C=  [0 ;  c2  ;  0  ;  0  ]  ; 
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D=[0;d2;0;d4]  ; 

K=  [0  0  0  0;  -15  0  30  0.3]; 

K=[0  0  0  0;  15  0  —30  —0.3];  %  corrected  by  Guojian  Lin 
%A=A+B*K; 

%cavity  model  information 

%osigma  = .  03*75  " 2/V " 2;  %sp  ecific  cavitator  data  from 

original  model 

%creates  input  file  for  scax  code 
cavlen=spline  (  sigmafit  ,lengthfit  ,  sigma); 

fid=fopen(  ’body  .  dat  5  ,  ’w’ )  ; 
fprintf(fid  ,  ’ body  .  dat \nCONE\nC\nD\n  ’ )  ; 
fprintf(fid  ,  ’32\tl6\t  — >MBOD,  _MEND\  n  ’ )  ; 

fprintf  (  fid  ,  ’  90 . 0 \  tO  .  5\ t->.HALF.(X>NE^\NGLE.  (DEGREES)  ,  .HEIGHT 
.OF  .ENDPLATE\  n  ’ )  ; 

fprintf(fid  ,  [  num2str  (  cavlen  )  ’ \t0 . 05\ t—>. CAVITY  .LENGTH,  . 

CAVITY  NODE,  FACTOR,,  (CJPT)  .TO.DEFINE.#.OF.NODES.IN.CAVITY\ 

n  ’  ] ) ; 

fprintf  (fid  ,  ’  l\t— >. ITERATE  .ON.CAVITATIONDTAIBER:  .1=ITERATE,  . 
0=NOT.ITERATE\  n  ’ )  ; 

fprintf(fid  ,[  num2str  ( sigma )  ’  \  t5 . 0E—5\t—>.CA VITATION .NUMBER. 
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TO.BE.OONVERGED.ON,OONVERGENCE.CRITERION\n  ’  ] )  ; 
fprintf  (  fid  ,  ’  1\ tO  .  5  \  t->.l— IF  .NON-DIMENSIONAL .WITH .BASE. 

DIAMETR.  .BASE.RADIUS.OTHERWISE\n  ’ )  ; 
fprintf(fid  ,  ’2 .  OE— 4\t— >^EPS^  (EPSILON)  ^fo  r  ^cavity,  length  \  n  ’ ) 
fprintf  (  fid  ,  ’  0 \ t  — > ^GRAVITY ^EFFECT :  ,1=AI)D,  ,0=IX),  NOT,  ADD\ii  ’ )  ; 
fprintf  (fid  ,  ’  104.0  \t— >J7ROUDE  JNfUMBER\n  ’)  ; 

fprintf  (fid  ,  ’  0\t->^GLEX>FAVITACK.EFFECT:  .1=ADD,  ,0=1)0,  NOT 
ADD\n  ’ )  ; 

fprintf  (  fid  ,  ’  1 0  \  t  — >^ANGLE^OF  OVTTACK^  ( IN  ^DEGREES)  \n  ’ )  ; 
fprintf  (fid  ,  ’  3 . 0E7\t— >.  REYNOLDS  JNUMBER’ )  ; 
fclose  (  fid  )  ; 

! scax . exe 

fid=fopen(  ’control  .  xls  ’)  ; 

Cd=str2num (fgetl(fid))  ; 

Mbod=str2num  ( fgetl(fid)); 

Mcav=str2num  ( fgetl(fid)); 
counter  =  1; 
while  ( 1 ) 

readstr=fgetl(fid)  ; 
if  readstr~=— 1 
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tester  (counter  ,:)  =str2num(  readstr  )  ; 
counter=counter  +  1; 

else 

break 

end 

end 

fclose  (  fid  )  ; 

%{ 

Cd  =  0. 85167; 

Mbod=32; 

Mcav=203; 
load  tester . txt 
%} 

xpos=tester ( : , 1 ) ; 
ypos=tester ( : ,2) ; 

xcp  =  [xpos  (Mbod+1) ;  tester  (Mbod+l:Mcav+Mbod, 3)  ]  ;  %adds  first 
x,y  value  to  the  cavity  section  of  xcp/ycp  data 
ycp  =  [ypos  (Mbod+1) ;  tester  (Mbod+l:Mcav+Mbod, 4)  ]; 

%xcp  =[x  ( Mbod+1) ;  xcp  ] ; 

%ycp  =[y  (Mbod+1) ;  yep]; 

normL— L/ (Rn*2)  ;  %finds  the  normalized  value 
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for  L 


index=max  ( find  (xcpCnormL)  )  ; 

xcptrunc  =  [xcp  ( 1 :  index  )  ;  normL  ]  ;  %integrates  only  up 

to  length  =  L 

ycptrunc  =  [ycp  ( 1 :  index  )  ;  interpl  (xcp  ,  yep  ,  normL)  ]  ; 

%  shift  C  o  eff=—Cd/8*  trap z  ( xcptrunc  ,  1 .  /  (  ycptrunc  .  " 2 )  )  ;  % 

normalized  shift  value  w/o  the  sin  ( alpha )*  cos  ( alpha )  term 
%unshiftedR=interp  1  ( xcp  ,  yep  ,  normL )  ; 

individualshift=— 1 Cd/8*  cumtrapz  (xcptrunc  ,  1 .  /  (ycptrunc. "2)  )  ; 

wmmmmmmmND  preamble 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


global  s  xO  xf  n  m  xf_index  obst_cords  obst_rs  lowest 
lowest  holder 
lowest  =  inf  ; 
lowest h o  1  d e r  =  [ ] ; 
s  =  14; 
m=2; 

x0  =  [0  0  0  0  0]’; 
xf  =  [— 20  0  0  0  80]’; 
xf _index  =  [1  5 ]  ; 
n=length  (xO )  ; 
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u.limit  =[pi  / 2 . 5 ]  ; 

obst_cords  =  [— 5  40]  \%coordinate  of  obstacle  (z,x) 
obst  _rs  =  [6 . 5] ; 

obst  _rs=obst  _rs  .  *  2  ;  %squares  the  radius  for  easier 
comparison  in  function 

guessT  =  1.1; 
guess.c  =0.6; 
guess.e  =0; 

guess  =  [ones(s  ,l)*guess_e  ones  ( s  ,l)*guess_c  ;  guessT  ]  ; 

%data  from  previous  run 

opt  ions=optimset  (  ’Algorithm  ’  ,  ’interior  —point  ’  ,  ’  MaxFunEvals  ’ 
,80000  ,  ’TolCon  ’  ,  le—  3)  ;%for  f  mine  on 
options=optimset  (  ’MaxFunEvals  ’  ,10000)  ;%for  fminsearch/ 
fminunc 

UB=[ones  (  s*m,  1 )  *  u_limit  ;  Inf  ]  ; 

LB=[— ones  (  s*m,  1 )  *  u_limit  ;  0  ]  ; 

%/or  lo op er  =  1 : 10 

tic 
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%[vars  fval]=  fminsearch  (  @  obj  ectiv  e_div  eplane  ,  guess  ),  options 


); 


Ainput  =  [eye  ( length  (  guess  ) ) ;  —  eye  ( length  (guess  )  )  ]  ; 
binput  =  [UB;—  LB]  ; 

for  looper=l:4 

[vars  fval]=  patternsearch  (  @obj ective.diveplane  ,  guess  , 
Ainput  ,  binput )  ; %,  options )  ; 
if  fval  <5 

guess=vars  ; 

obst_rs=(sqrt(obst_rs)+.5)  "2; 

else 

break 

end 

end 

toe 

%guess=vars ; 

%end 
obst  _rs 
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for  i  =  1  :m 


u  ( :  ,  i  )=vars  ( ( i  —  l)*s+l:s*i)  ; 

end 

control=u ; 

T=vars  (end)  ; 

control 

T 


delta_e  =  [];  delta.c  =  []; 

[  i  n  t  _  t  ,  int_x]  =  ode45  (  @dive_plane_noncyl  ,  [0  T/s],x0,[], 
control  (1  ,:)  )  ; 

delta_e=[delta_e  ;  ones  ( length  ( in  t  _t  )  ,l)*control(i  ,1)  ]; 
delta_c=[delta_c  ;  ones  ( length  ( in  t  _t  )  )  *  control  ( i  ,2)  ]  ; 

times=T/ s  :T/s  :T; 
for  i  =2:  length  ( times  ) 

[t,  x]  =  ode45  (  @dive_plane_noncyl  ,[  times  ( i  — 1)  times(i)], 
i  n  t  _x  ( end  ,:)  ’,[],control(i,:)); 
int  _t  —  [  int  _t  ;  t  ]  ; 
int _x  =  [  int _x  ;  x  ]  ; 

dclta_e=[delta_e  ;  ones  ( length  (t)  ,l)*control(i  ,1)  ]; 
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delta_c=[delta_c  ;  ones  ( length  (t)  ,l)*control(i  ,2)  ]; 

end 

K=[0  0  0  0;  0  0  0  —0.9];  %inner  loop  control  law 
tot_control=([delta_e  delt  a_c  ’]  +K*  int  _x  ( :  ,  [  1  2  3  4] )  ’ )  ’ ;  % 
makes  columns  with  [delta-e-tot  ,  delta-C-tot] 
delta_c_tot=tot_control  (:  ,2)  ; 


figure  ( 1 ) 

plot ( int  _x ( :  , 5 )  ,  int  _x  ( :  , 1 ) ) ; xlabel ( ’x^  (m)  ’ ) ; ylabel ( ’ z „  (m)  ’ )  ; 

hold  on 

[t.inneronly  ,  x.inneronly]  =  ode45  ( @dive_plane_noncyl  ,  [0  T]  , 
xO  ,  []  ,  zeros  (m, 2)  )  ; 

plot(x_inneronly(:  ,5)  ,x_inneronly(:  ,1)  ,  ’ — b  ’ )  ;  legend  (  ’  inner  ^ 
and^outer  ^loop  ’  ,  ’  inner  ^loop  ^only  ’ ) 

%plots  the  obstacles 
for  i  =1:  length  (  obst  _rs  ) 
xc=obst _cor ds ( i  ,2)  ; 
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zc=obst _cor ds  ( i  ,2)  ; 
rad=sqrt  (obst_rs(i)); 
xes=linspace(— rad  ,  rad)  ; 

z _plus=sqrt  ( rad "2— xes  .  ~  2  )+obst  _cords  ( i  ,1)  ; 
z_minus=—  sqrt  (  obst  _rs  (  i  )— xes  .  '  2  )+obst  .cords  ( i  ,  1 )  ; 
plot ( xes+xc , z.plus  ,  ’ r ’ )  ; plot ( xes+xc , z .minus ,  ’ r  ’ ) ; 

end 

legend  (  ’inner  ^and^  outer  ^loop  ’  ,  ’inner^loop^only  ’  ,  ’obstacle  ’) 

%pl  ottin  g  the  body  orientations 

indexes  =  [find(int_t>int_t  (end)  *.20,1)  find(int_t>int_t  (end) 
*.40,1)  find  ( int  _t  >int  _t  (end)  *  .  65 , 1 )  length  (  int  _t  )]  ; 
ef  =  3 -^enlargement  factor 
I^L*  ef ;R=R*  ef ; 

xcordbody  =  [linspace  (0 ,  —  L  , 4 )  linspace(— L ,  0  ,4)  ]  ; 
zcordbody  =  [0  R  R  R  — R  — R  — R  0]  ; 

for  i  =1:  length  ( indexes  ) 

theta=int_x  ( indexes  ( i )  ,3)  ; 
xpos=int  _x ( indexes ( i )  ,5)  ; 
zpos=i nt  _x ( indexes ( i )  ,  1 ) ; 
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x  cor  d=xpos+xcordbody  *  cos  (  —  theta  )— z  cor  dbody*sin(—  theta)  ; 
zcord=zpos+zcordbody  *  cos  (  —  theta  )— (— xcordbody  )*sin(—  theta 


patch  (xcord,zcord,  ’  r  ’ ) 

w  =  int_x(indexes(i)  ,2)  ; 

alphacav=w/V+delta_c_tot  ( indexes  ( i ) )  ; 

cavitycordstop  =  (ycptrunc+individualshift*sin(  alphacav)* 
cos  (  alphacav  )  )  *2*Rn; 

cavitycordsbottom=(—  ycptrunc+in  divi  duals  hift*sin( 
alphacav  )  *cos  (  alphacav  )  )  *2*Rn; 

xcords=— xcptrunc *2*Rn; 

cavitycordstop=cavitycordstop  *  ef  ;  cavitycordsbottom= 
cavitycordsbottom*ef  ;  xcords=xcords *  ef  ; 

cavxt=xpos+xcords  *  cos  (  — theta+w/V)  — cavity  cor  dst  op  *  sin  (— 
theta+w/V)  ; 

cavzt=zpos+cavitycor  dst  op  *cos(—  thet  a+w/V)  —  (—  xcords)*sin 
(—theta+w/V)  ; 

cavxb=xpos+xcords  *  cos  (  —  thet  a+w/V)  — cavity  cor  dsbot  tom  *  sin 
(—theta+w/V)  ; 


165 


cavzb=zpos+cavitycordsbottom  *cos(—  t  het  a+w/V)  —  (—  xcords  )  * 


sin(— theta+w/V)  ; 

plot ( cavxt ,cavzt  ,  ’  g  ’ ) ;  plot ( cavxb , cavzb , ’ g  ’ ) 

end 

L=L/ef ;R=R/ef ; 


V/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/OA 


figure  ( 2 )  ; 

subplot  (2,2,1)  ;plot(int_t  ,  int_x  (:  ,1)  )  ;  y  label  (  ’  depth  ^z  ^  (m)  ’ )  ; 
grid  on  ;  xlabel  (  ’  time  (  s  )  ’ )  ; 

subplot  (2,2,2)  ;  plot  ( int  _t  ,  int  _x  (:  ,2)  )  ;  y  label  (  ’velocity  (m/ 
s  )  ’ )  ;  grid  on  ;  xlabel  (  ’  time  ^  (  s  )  ’ )  ; 

subplot  (2,2,3)  ;  plot  ( int  _t  ,  int  _x  (:  ,3)  )  ;  y  label  (  ’  pit  elm  angle 
theta^(rad)  ’)  ;grid  on;  xlabel  (  ’  time  ^  (  s  )  ’ )  ; 

subplot  (2,2,4)  ;  plot  ( int  _t  ,  int  _x  (:  ,4)  )  ;ylabel(  ’pitch  ^rate^q^( 
rad /  s  )  ’ )  ;  grid  on  ;  xlabel  (  ’  time  u  (  s  )  ’ )  ; 

figure  ( 3 )  ; 

subplot  (2 ,1  ,1)  ;  plot(int_t  ,  de It a_e  )  ;ylabel(  ’Fiimangle^rad)  ’) 

;  title  (  ’  Outer  ^  loop  ^  control  ^only  ’ )  ;xlabel(  ’  time  ^  (  s  )  ’)  ; 
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subplot  (2 ,1  ,2)  ;  plot(int_t  ,  de It a_c  )  ;ylabel(  ’  Cavitator^angle^( 
rad)  ’ )  ;  title  (  ’Outer^loop^control^only  ’ )  ;  x  lab  el  (  ’  time  ^  (  s  )  ’ 


figure  (4) 

subplot  (2,1,1)  ;plot(int_t  ,delta_c_tot);ylabel(’Cavitator^ 
angle  „  (  rad )  ,  wtot  al  ’ )  ;  xlabel  (  ’  1  i  me  s  )  ’ )  ; 
subplot  (2,1,2)  ;  plot  ( in  t  _t  (2:  end)  ,diff(delta_c_tot)./  diff( 
int  _t ) )  ;ylabel(  ’Cavitator^angle^rate^(rad/s)  ,  ^  t  ot  al  ’)  ; 
xlabel  (  ’  time  ^  (  s  )  ’ )  ; 

function  val  =  obj  ect  i ve _di veplane  ( var  ) 

%var  is  split  up  into  [xl(ts)  x2(ts)  ...  xn(ts)  ul(ts)  u2(ts 
)  ...  um(  ts )  T] 

%the  ts  are  split  up  to  T/s,  s  is  the  number  of  steps 
%n  is  the  number  of  states  ,  m  is  the  number  of  control 
v ariab  les 

%ox0  are  the  initial  states,  xf  is  the  final  state  and 
%oxf -index  is  the  index  of  the  final  states  that  are  desired 
to  be  fixed  . 

%for  instance  if  only  the  first  two  states  are  desired  to  be 
fixed  , 
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%xf -index  =[1  2] 


global  s  xO  xf  m  xf.index  obst.cords  obst_rs  lowest 
lowest  holder 

for  i  =  1  :m 

u  ( :  ,  i  )=var  ( ( i  —  l)*s  +  l:s*i)  ; 

end 

%  creates  a  mx  of  x  =  [ul(l)  u2  ( 1 ) 

%  ul (2)  u2  (2) 

% 

%  ul ( s )  x2  ( s ) 

T=var (end)  ; 

c  =  []; 

totalx  =  [];  %keep  track  of  entire  trajectory 

%equality  constraints :  dynamic  constraints  betweens  states, 
and  gluing 

%constraint  for  x(s)=xf 


un(  1)  ; 
un(2)  ; 

un( s ) ; 
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%constraint  for  first  x(l) 

[  t  xtemp]  =  ode45  (@dive_plane_noncyl  ,  [  0  ,  T/s  ]  ,  xO  ,  [  ]  ,  u  ( 1  , : )  )  ; 
x2=xtemp  (end  ,  :  )  ; 
totalx  =  [totalx  ;  xtemp  ]  ; 

%constraints  to  make  sure  w<6 

for  i  =  2 :  s 

[  t  xtemp]  =  ode45  (@dive_plane_noncyl  ,  [T/s*(i— l),T/s*(i)]  , 
x2  ’  ,  [  ]  ,  u  ( i  , : )  )  ; 
x2=xtemp  (end  , :  )  ; 
c  =  [c;  abs  (x2  (2)  )  — 6] ; 
totalx  =  [totalx  ;  xtemp  ]  ; 

end 

%ob  stacle  constraints 
for  i  =1:  length  (  obst  _rs  ) 

closest  d  is  t=min(  ( t  ot  alx  ( :  ,1)  — obst_cords(i  ,  1 ) )  .  ~  2+ ( ( 
t  ot  alx  ( :  ,5)  —  obst_cords(i  ,  2 )  )  .  "  2 )  )  ; 
c  =  [c;  obst.rs  ( i )  — closestdist  ]  ; 

end 

c 
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%control  bounds 


for  i  =1:  length  (  var )  — 1 

c  =  [c;  abs  (  var  ( i ) )— pi /2 . 5] ; 

end 

%enforces  final  condition 
ceq=xf  (  xf.index  )— x2  (  xf.index  )  ’ ; 
u 

T 

c 

ceq 

val=var(end)  +  100*sum(  (  c  >0)  .  *  abs  (  c  )  )  +100*sum(  abs  (  ceq )  ) 
val 

i f  valclowest 
lowest=val  ; 
lowest  hoi der=var  ; 

end 

function  dxdt  =  di  ve_plane_noncy  1  ( t  ,  y  ,  u) 

%  f— state  model  from  Dzielski  &  Kurdila  ’s  paper 
%  x=[z  w  theta  q] 
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global  Rcdot  R  V  L  A  B  C  D  Rn  Re  gcounter  individualshift 


ycptrunc  xcptrunc ; 

x=y  (1:4,:); 
w  =  x(2  ,:)  ; 
delta_e  =  0; 

%delta_c  =  —  15*x  ( 1  , : )  +30*x  ( 3  , : )  +0.3*x  (4  ,  •' )  ;  (Kurdila  in  JVC) 
delta_c  =  15*x  ( 1  , : )  —  30*x  ( 3  , : )  fe0.3*x  (4  , : )  ;  %—2*x  (2  ,  :)/V;  % 
correction  by  Guojian  Lin 
delta_c  =— 0.9*x  (4  ,:)  ; 

%delt a-C  =30*x  ( 1  ,:)  —  60*x  ( 3  ,:)  —0.6*x  ( 4  ,'■);% applying  true  2x 
feedback  ( affects 

%cavity  shape  as  well  instead  of  just  adding  it  into  the  A 
mx  with 

%A=A=B*K  (where  K  contains  the  linear  feedback  law) 

delta_c  =  delt  a_c+u  (2 )  ; 
delta.e  =  delt  a_e+u  ( 1 )  ; 

%delta_c  =  —30*x(3,:)—0.3*x(4,:);%—2*x(2,:)/V;  %no  depth 


171 


feedback 


%od  elt  a  _c  =0  %no  feedback; 

%odelta  -C  =15*x  ( 1  , : )  —  30*x  (3  , : )  —  6*x  (4  ;%b  ifurcatio  n  control 

%odelta_c  =  1 5*  x  ( 1  , : )  —90*x  (  3  , : )  —  0. 9*  x  ( f  , : )  ;  %gain  augmented 
feedback 

alphacav=w/V+delt a_c  ; 


%  if  (w>0)  %planes  on  bottom  Rcdot<0  so  sgn  ( alpha )—sgn  (w) 

%  alpha  =  w/V; 

%  else  %planes  on  top 
%  alpha  =  w/V; 

%  end 

alpha=atan  (w/V)  ; 

if  (w>0) 

cavity  cor  ds=abs(— ycptrunc+indi  vi  du  al  s  hi  ft  *sin  (  alphacav  )  * 
cos  (  alphacav  ) )  *2*Rn;  %unnormalized  cavity  radii 
along  planing  location 

else 
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cavitycords=abs  (  ycptrunc+ individuals  hi  ft  *sin(alphacav)* 
cos  (  alphacav  ) )  *  2  * Rn ; 

end 

delt  a=abs  (  ycptrunc  *2*Rn)— R;  %individual  deltas  ( 

difference  in  radius  of  cavity  vs  body)  unnormalized 
xcords=xcptrunc  *2*Rn ;  %unnormalized  x  p  ositions 

%alphag  eometric  =  atan(w/V); 
body cor ds=abs (xcords*t an (alpha) )+R; 

%r efining  mesh  only  where  planing  begins 

index=find  ( (  bodycords>cavitycords  )  .  *  (  xcords>max(  xcords)  /  3) 
,1);  %finds  first  location  of  planing  of  the  cylinder 
part  of  the  body 
if  index>0  %if  planing 

xcordsnew=linspace  (  xcords  ( index  —2)  ,  xcords  (end)  )  ’ ; 
cavitycords=interpl  (xcords  ,  cavitycords  , xcords  new)  ; 
delt  a=abs  ( interpl  ( xcords  ,  ycptrunc  ,  xcordsnew )  *2*Rn— R)  ; 
xcords=xcordsnew  ;  %xcords  becomes 

truncated  to  just  the  planing  area 
body cor ds=abs (xcords*t an (alpha) )+R; 

end 
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hdepths=(  body  cords—  cavity  cords  )  .  *  (  bodycords>cavitycords  )  ; 

%finds  immersion  depths  ( only  where  body  is  planing 

) 

if  index>0  %only  does  this  calculation  if  planing 

for  j  =1:  length  ( xcords  )  %requires  ’’for”  loop  since 

delta+hdepths  can  =0  giving  NaN 
if  ( bodycords  ( j  )>cavitycords  ( j  ) )  %calculation 
only  where  planing  (on  now  truncated  xcords) 
integrand  (j  )  =2*  delta  (j)"2./(delta(j  )+hdepths  ( j  ) ) 
3 ; 

else 

integrand  ( j  )  =0; 

end 

end 

else 

integrand  ( 1 :  length  (xcords)  )=0; 

end 


hdots=V*sin  ( abs  (alpha )— atan  (diff(cavitycords)  ./  diff(  xcords)) 
)\%cavity  angle  =  atan  (  diff  (  cavity  cor  ds  )./ diff  (  xcords  )  ) 
%ohdots=V*  sin  ( alpha )  ; 
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FpnormNoncyl=— sum(  diff  (cumtrapz  (xcords  ,  integrand  ’)  )  .*hdots 
. '  2 )  *sign  (  alpha  )  ; 

%this  DOES  take  into  account  the  cavity  slope,  d  if  f  ( cumtrapz 
)  gives  the 

%individual  trapezoidal  areas,  this  multiplied  by  hdot/'2 
over  that 

%particular  area  gives  the  individual  integral  totals, 
summed  up  gives  the 

%entire  planing  force. 

%FpnormNoncyl——V" 2*  trap z  ( xcords  ,  integrand )*  sin  ( alpha )  "2*  sign 
(alpha);  %sign  to  get  the  proper  sign  on  planing  force 

%this  also  assumes  that  the  planing  angle  is  alpha  ( does  not 
take  into 

%account  slope  of  cavity  shape).  also  alpha  in  this  case 
takes  into 

%account  the  Rcdot  term,  is  that  right?? 

if  FpnormNoncyl~=0 

%xp=—  V*2.  *  trap  z  ( xcords  ,  integrand  .*  xcords  )  .*  sin  (  alpha  ) 

* 2 .*  sign  ( alpha )./ FpnormNoncyl ;  %effective  point 
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location  of  planing  force  (from  front  of  vehicle) 
xp=— sum(  diff  ( cumtrapz  ( xcords  ,  integrand  ’ .  *  xcords  ) )  .  *  hdots 


.  ~  2  *  sign  (  alpha )  )  .  /  FpnormNoncyl ; 

else 

xp=L; 

end 


gcounter=gcounter  +  1; 


n=FpnormNoncyl ; 


dxdt_4  =  A*x+B*  [  dclta_e  ;  delta_c]+C-|-D*n  .  *  [1  1  1  (xp)/L]’; 


theta=x  ( 3  , : )  ; 

dzdt^w*cos  (theta  )— V*sin  (theta  )  ; 
xdot=V*cos (theta )+w*sin  (theta ) ; 

dxdt  =  [dzdt;  dxdt_4([2  3  4]  ;  xdot  ]  ; 

B.3  Maneuvering  with  Delay  and  Impact  Model 
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clc 


clear 
close  all 


global  Rcdot  RpRVLABCD  delay  V; 

mmmmmmmmmmmmc onst ants  section 


g  =  9.81; 
m  =  2 ; 

Rn  =  0.0191; 

R  =  0.0508; 

L  =  1.8; 

sigma  =  0.0241;  %0.035  or  0.03 

sigma  =  0. 03; 
n=  0.5; 

CxO  =  0.82; 

Cx  =  CxO*(l  +  sigma)  ; 

C  =  l/2*Cx*RiA2/R7  2; 

V  =  sqrt  (0. 03*75  ~  2/ sigma)  ; 

M0=  [7/9  17*L / 3 6 ;  17*L/36  1 1/60*R7 2  +  133/405*L  '  2 ] ; 

A0=C*V*[(1  —  n )  /m/L  —n/m;  — n/m  — n*L/m]+V*  [0  7/9;  0  17*L/36] 
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B0=C*V~2*[  —  n/m/L  1/m/L;  — n/m  0]; 

C0=  [7/9 ;  17*L /36]  * g  ; 

D0=  [1/m/L;  1/m]  ; 

D00  =  [l/m/L  0;  0  1/m/L]; 

%  correction  by  Guojian  Lin 

A0=C*V*[(  —  1  —  n)/m/L  —n/m;  —  n/m  —  n*L/m]+V*  [0  7/9;  0  17*L/36]; 
B0=C*V~2*[  — n/m/L  —1/m/L;  —n/m  0]; 

A2=inv  (MO)  *A0 ; 

B2=inv  (MO)  *B0 ; 

C2=inv  (MO)  *C0 ; 

D2=inv  (MO)  *D0 ; 

D22=inv  (MO)  *D00 ; 

K1  =  L/Rn/  ( 1 . 92/  sigma  —  3)  —  1; 

K2  =  sqrt  (1  —  (1  —  4.5* sigma/(l  +  sigma )  )  *K1 "  (40/ 1 7)  )  ; 

Rc  =  Rn* (0. 82* ( 1  +  sigma)  /sigma)  "0.5* K2; 

Rp  =  (Rc— R)  /R; 

Rcdot  =  —20  /17*  (0.82*  (1  +  sigma )  /  sigma )  '  0 . 5*V*  (1  —4.5*  sigma  /(In¬ 
signia  ) )  *K1  (23/17)  /K2 / ( 1 . 92/ sigma  —3) ; 
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a22=A2  (1,1);  a24=A2(l,2);  a42=A2(2,l);  a44=A2(2,2); 
b21=B2 (1,1);  b22=B2 (1,2);  b41=B2(2,l);  b42=B2(2,2); 
c2=C2(l,l); 

d2=D2 (1,1);  d4=D2(2,l); 

d21=D22  (1,1)  J  d22=D22(l  ,2)  ;  d41=D22  (2  , 1 )  ;  d42=D22  (2  ,2)  ; 

A=[0  1  -V  0;  0  a22  0  a24 ;  0  0  0  1;  0  a42  0  a44  ]  ; 

B=  [0  0;  b21  b22 ;  0  0;  b41  b42  ]  ; 

C=  [0 ;  c2  ;  0  ;  0  ]  ; 

D=[0;d2;0;d4]  ; 

K=  [0  0  0  0;  -15  0  30  0.3]; 

K=[0  0  0  0;  15  0  —30  —0.3];  %  corrected  by  Guojian  Lin 
K=[0  0  0  0;  0  0  0  —0.9];  %changed  to  remove  dependence  on  z 
or  theta 

%A=A+B*K;  %this  has  the  built  in  inner  loop  control  already 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


delay=L/V ; 


global  s  xO  xf  num.in  xf.index  obst_cords  obst.rs  u_init 
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s  =  14; 


num  _in  =  2; 
x0  =  [0  0  0  0  0]’; 
xf  =  [— 20  0  0  0  80]’; 
xf .index  =  [1  5 ]  ; 

%n=length  (xO ) ; 
u_limit  =[pi  / 2 . 5 ] ; 

obst_cords  =  [— 5  40]  \%coordinate  of  obstacle  (z,x) 
obst_rs  =  [6]; 
u  _i  ni  t=zeros  ( num_in  ,  1 )  ; 

obst  _r s=obst  _r s  . ~  2  ;  %squares  the  radius  for  easier 
comparison  in  function 


guessT  = .  8 ; 
guess.c  =0; 
guess.e  =0; 

guess  =  [ones(s  ,l)*guess_e  ;  ones  ( s  ,l)*guess_c  ;  guessT  ] 

options  =  psoptimset  (  ’TolCon  ’  ,le  —  3)  ; 

UB=[ones  (  s*num_in  ,  1 )  *  n.limit  ;  Inf  ]  ; 
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LB=  [ —  ones  (  s*num_in  ,  1 )  *  u_limit  ;  0 . 1  ]  ; 

tic 

%[vars  fval]=  fmincon  (  @  obj  ectiv  e_diveplane  ,  guess 
>[]>[]>[]>[]  .  LB,  UB,  [  ]  ,  options)  ; 

[vars  fval  exitflag]=  patternsearch  (  @obj  ective_diveplane  , 
guess  ,  []  ,  []  ,  []  ,  []  ,  LB ,  UB ,  [  ]  ,  options)  ; 

toe 


for  i=l:num_in 

u  ( :  ,  i  )=vars  ( ( i  —  l)*s+l:s*i)  ; 

end 

control=u ; 

T=vars  (end)  ; 

control 

T 

for  i=l:num_in 

upp(  i  )=spline  (0  :T/s  :T,  [  u_init  (  i  )  ;  u  ( :  ,  i  )  ] )  ; 

end 
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sol  =  dde23 (@( t , y , Z)  dive_plane (t , y , Z , upp)  ,[ delay ]  ,@(  t )  [xO 


(1:4)  ;  sqrt  (V"2+x0  (2)  "  2)  *t  ]  ,[0  T]  ,[])  ; 
i  n  t  _  t  =  s  o  1  .x  ’  ; 
int  _x=s o  1  .  y  ’  ; 


figure  ( 1 ) 

plot ( int  _x ( :  , 5 )  ,  int  _x  ( :  , 1 ) ) ; xlabel ( ’x^  (m)  ’ ) ; ylabel (  ’  z , ,  (m)  ’ ) 

hold  on 

u_inner=zeros  ( s  ,m)  ; 
for  i=l:num_in 

upp.inner  (  i  )=spline  ( 0  :T/s  : T ,  [  u  _in i t  (  i  )  ;  udnner  ( :  ,  i  )  ] )  ; 

end 

sol.inner  =  dde23  (@( t  ,  y  ,  Z)  dive.plane  ( t  , y  ,  Z  ,  upp.inner  )  ,[ 
delay  ]  ,@(  t )  [xO  ( 1 : 4)  ;  sqrt  (\T2+xO  (2)~2)*t],[0  T],[]); 
t_inneronly=sol_inner  . x  ’  ; 
x_inneronly=sol_inner  .y’; 

plot  (  x  _  i  n  n  e  r  o  n  1  y  ( :  ,  5  )  ,  x  _  i  n  n  e  r  o  n  1  y  ( :  ,  1 )  ,  ’  b —  ’ )  ; 
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%plots  the  obstacles 
for  i  =1:  length  (  obst  _rs  ) 
xc=obst _cords  ( i  ,2)  ; 
zc=obst _cords  ( i  ,2)  ; 
rad=sqrt  (obst_rs(i)); 
xes=linspace(— rad  ,  rad)  ; 

z_plus=sqrt(rad~2— xes  .  ~  2  )+obst  _cords  ( i  ,  1 )  ; 
z_minus=—  sqrt  (  obst  _rs  (  i  )— xes  .  '  2  )+obst  _cords  ( i  ,  1 )  ; 
plot  ( xes+xc  ,  z_plus  ,  ’ r ’ )  ;  plot  (  xes+xc  ,  z_minus  ,  ’  r  ’ )  ; 

end 

legend  (  ’inner  ^and^  outer  ^loop  ’  ,  ’inner^loop^only  ’  ,  ’obstacle  ’) 

figure  ( 2 )  ; 

subplot  (2,2,1)  ; plot ( int  _t  ,  int  _x  (:  ,1)  )  ;  y  label  (  ’  depttmz  ^  (m)  ’ )  ; 
grid  on  ;  xlabel  (  ’  time  ^  (  s  )  ’ )  ; 

subplot  (2,2,2)  ;  plot  ( int  _t  ,  int_x  (:  ,2)  )  ;  y  label  (  ’velocity  (m/ 
s  )  ’ )  ;  grid  on  ;  xlabel  (  ’  time  ^  (  s  )  ’ )  ; 
subplot  (2,2,3)  ;  plot  ( int  _t  ,  int_x  (:  ,3)  )  ;  y  label  (  ’  pit  elm  angle 
theta  ^  (  rad )  ’ )  ;  grid  on  ;  xlabel  (  ’  time  ^  (  s  )  ’ )  ; 
subplot  (2 ,2  ,4)  ;  plot  ( int  _t  ,  int  _x  (:  ,4)  )  ;  y  label  (  ’pitch  urateuqu( 
rad /s  )  ’ )  ;  grid  on  ;  xlabel  (  ’  time  u  (  s  )  ’ )  ; 
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figure  ( 3 )  ; 

delta_c  =  ppval  (upp  (2)  ,  int  _t  )  ; 
delta.e  =  ppval  (upp  ( 1 ),  int  _t )  ; 

subplot  (2 ,1  ,1)  ;  plot  ( int  _t  ,  de It a_e  )  ;ylabel(  ’Fiimangle^(rad)  ’ ) 
;  title  (  ’  Outer  ^  loop  ^  control  ^only  ’ )  ;xlabel(  ’  time  ^  (  s  )  ’ )  ; 
subplot  (2 ,1  ,2)  ;  plot ( int  _t  ,  de  It  a_c  )  ;ylabel(  ’  Cavitatoruangle^( 
rad)  ’ )  ;  title  (  ’Outer^loop^control^only  ’ )  ;  x  lab  el  (  ’  time  ^  (  s  )  ’ 


figure  (4) 

tot_control=([delta_e  ’ ;  delt  a_c  ’  ]  +K*  int  _x  ( :  ,  [  1  2  3  4] )  ’ )  ’ ;  % 
makes  columns  with  [delta-e-tot  ,  delta-C-tot] 
delta_c_tot=tot_control  (:  ,2)  ; 

subplot  ( 2 , 1  , 1 )  |;plot  ( int  _t  ,delta_c_tot);ylabel(’Cavitator^ 
angle  ^  (  rad )  ,  ^  tot  al  ’ )  ;xlabel(  ’  time  ^  (  s  )  ’ )  ; 
subplot  (2,1)2)  ;plot(int_t  (2:  end)  ,diff(delta_c_tot)./  diff( 
int  _t ) )  ;ylabel(  ’Cavitator^angle^rate^(rad/s)  ,  ^-t  ot  al  ’ )  ; 
xlabel  (  ’  time  ^  (  s  )  ’ )  ; 

function  val  =  obj  ect  i  ve_di  veplane  ( var  ) 

%var  is  split  up  into  [ul(ts)  u2(ts)  ...  um(ts)  T] 

%xf -index  is  the  index  of  the  final  states  that  are  desired 
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to  be  fixed  . 

%for  instance  if  only  the  first  two  states  are  desired  to  be 
fixed  , 

%xf  -index  =[1  2] 

global  s  xO  xf  num.in  xf.index  obst_cords  obst.rs  u_init 
delay  V 

for  i=l:  num.in 

u  ( :  ,  i  )=var  ( ( i  —  l)*s  +  l:s*i)  ; 

end 

%  creates  a  mx  of  x  =  [ul(l)  u2 ( 1 ) 

%  ul  (2)  u2(2) 

% 

%  ul(s)  x2(s) 

T=var (end) ; 

c  =  []; 

totalx  =  [];  %keep  track  of  entire  trajectory 

for  i=l:  num.in 


un(  1)  ; 
un(2) ; 

un( s )  ; 
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upp  (  i  )=spline  ( 0  :T/s  :  T ,  [  u  _i  ni  t  (  i  )  ;  u  ( :  ,  i  )  ] )  ; 

end 

%equality  constraints :  dynamic  constraints  betweens  states, 
and  gluing 

%constraint  for  x(s)=xf 

%constraint  for  first  x(l) 

sol  =  dde23  (@(  t  ,  y  ,  Z)  dive_plane ( t , y , Z , upp)  ,[ delay ]  ,@(  t )  [xO 
(1:4)  ;  sqrt  (V~2+x0  (2)  *2) *t ]  ,  [0  T]  ,[])  ; 
t=s  o 1 . x ; 
totalx=sol  .  y  ’ ; 

%constraint  to  make  sure  w<6 
c  =  [c;  max(  abs  ( tot  alx  ( :  ,  2 )  )  )  — 6] ; 

%ob  stacle  constraints 
for  i  =1:  length  (  obst  _rs  ) 

closest  d  is  t=min(  ( t  ot  alx  ( :  ,1)  — obst_cords(i  ,  1 ) )  .  ~  2+ ( ( 
tot  alx  ( :  ,5)  — obst_cords(i  ,  2 )  )  .  "  2 )  )  ; 
c  =  [c;  obst_rs  ( i )  —  closestdist  ]  ; 

end 
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xend=tot alx (end , :  )  ; 

%enforces  final  condition 
ceq=xf  (  xf  _index  )— xend  (  xf  .index  ) 

val=var(end)  +  100*sum(  (  c  >0)  .  *  abs  (  c  )  )  +100*sum( abs  (  ceq )  ) 
u 

T 

c 

ceq 

val 

function  dxdt  =  dive  .plane  ( t  ,  y  ,  Z  ,  npp ) 

%  f— state  model  from  Dzielski  &  Kurdila  ’s  paper 
%  x=[z  w  theta  q  x] 

global  Rcdot  Rp  R  V  L  A  B  C  D; 

thetat=Z (3,1)  ; 

xt=Z (5,1)  ; 

zt=Z ( 1  ,1)  ; 

wt=Z (2,1)  ; 

thetact=thetat— atan(wt/V)  ; 
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z=y ( i , ; ) ; 

(  2  , : )  ; 

theta=y ( 3  , : )  ; 

q=y ( 4  , : )  ; 
x=y ( 5  , : )  ; 

oldstates=y(l:4  ,:)  ; 

zdot^w*cos (theta )— V*sin (theta )  ; 
xdot=V*cos (theta )+w*sin (theta ) ; 

a=cos ( thet act ) * (x— xt )+sin  ( the tact ) * ( zt— z ) ; 

a=L*cos ( theta— thetact ) ;  %this  is  necessary  when  dropping  the 
delay  close  to  0 

b=sin  (thetact)  *  (x— xt )— cos  (thetact)  *(zt— z)  ; 
adot=xdot *cos (thetact)— zdot*sin( thetact)  ; 
bdot=xdot *  sin  ( thetact )+zdot*cos ( thetact ) ; 

hhat=a*tan  ( theta— thet  act  )+b  ;  %this  is  the  relative  position 
of  the  body  wrt  cavity  ( can  be  neg  or  pos) 
h=max(  [  abs  (  hhat )—  Rp*R,  0  ] )  ;  %this  is  the  immersion  depth, 
positive  no  matter  what  side  is  immersed 
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hdot  =  (a*sec(theta— theta  ct)  "2*q+adot  *t  an  (theta—  thetact  )+bdot ) 
*sign  ( hhat )— Rcdot ;  %the  multiplication  by  sign(hhat)  is 
to  make  sure  that  the  immesion  rate  is  in  the  direction 
of  the  immersion 

hdot=max( hdot  , 0 )  ;  %the  planing  force  formulation  only  works 
for  an  immerstion  rate  that  is  in  the  direction  of  the 
immersion  so  only  positive  immersion  rates  are  available 


delta=Rp*R; 

i  f  h>0 

Fplane=—  hdot  "2*1/  tan  ( abs  (theta— thetact ) ) .  *  (1  —  delta  "2./  (h 
+  delta)  ,"2)*sign(hhat)  ; 

else 

Fplane  =0; 

end 

delta.c  =  ppval  (upp  (2)  ,  t )  —  0.9*  oldstates  (4  , : )  ; 
delta_e  =  ppval  (upp  ( 1 ),  t )  ; 
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n=Fplane  ; 


dolddt  =  A*  oldst  at  es+B*  [  delt  a_e  ;  delt  a  c ]+C+D*n  ; 


dxdt  =  [zdot  ;  dolddt  ( [2  3  4],:);  xdot  ]  ; 
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